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cn Abstract 

CN We calculate the leptonic decay constants of and D^^) mesons in lattice QCD using staggered 

I light quarks and Fermilab bottom and charm quarks. We compute the heavy- light meson corre- 

T-H lation functions on the MILC asqtad-improved staggered gauge configurations which include the 

J> effects of three light dynamical sea quarks. We simulate with several values of the light valence- and 

^ sea-quark masses (down to ~ ms/10) and at three lattice spacings (a « 0.15, 0.12, and 0.09 fm) and 

extrapolate to the physical up and down quark masses and the continuum using expressions derived 
in heavy-light meson staggered chiral perturbation theory. We renormalize the heavy-light axial 
current using a mostly nonperturbative method such that only a small correction to unity must be 
computed in lattice perturbation theory and higher-order terms are expected to be small. We obtain 
fB+ = 196.9(8.9) MeV, /b, = 242.0(9.5) MeV, /b+ = 218.9(11.3) MeV, /d, = 260.1(10.8) MeV, 
and the SU(3) flavor-breaking ratios /bV/b = 1.229(26) and /d,//d = 1.188(25), where the 
numbers in parentheses are the total statistical and systematic uncertainties added in quadrature. 

PACS numbers: 12.38.Gc, 13.20.Fc, 13.20.He 
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I. INTRODUCTION 



Leptonic decays of B and D mesons, in which the hadron annihilates weakly to aW boson, 
are important probes of heavy-to-light quark flavor- changing interactions. When combined 
with a nonperturbative lattice QCD calculation of the heavy-light pseudoscalar meson decay 
constant, fs or a precise experimental measurement of the leptonic decay width allows 
the determination of the Cabibbo-Kobayashi-Maskawa (CKM) quark-mixing matrix element 
\Vub\ or \Vcd\- Conversely, if the relevant CKM matrix element is known from an independent 
process such as semileptonic decay or from CKM-unitarity constraints, a comparison of the 
decay constant from lattice QCD simulations with that measured by experiment provides a 
straightforward test of the Standard Model. As the lattice and experimental determinations 
become more precise, this test will become more sensitive and may ultimately reveal, through 
the appearance of a discrepancy, the presence of new physics in the quark flavor sector. 

Improved determinations of the B meson decay constant fs are of particular importance 
given the current, approximately 3-cr tension in the CKM unitarity triangle that may indi- 
cate the presence of new physics in iJ^-mixing or S —> rz/ decay [1-3]. The experimental 
uncertainty in the branching fraction B{B — )■ tu) is at present ~ 30% [4, 5], but this error is 
expected to be reduced to ~ 10% by Belle II at KEK-_B in as little as five or six years [6, 7], 
at which point even modest improvements in the determination of Jb will significantly help 
constrain the apex of the CKM unitarity triangle and isolate the source of new physics [8]. 

Because leptonic decays are "gold-plated" processes in numerical lattice QCD simulations 
(they have a single stable hadron in the initial state and no hadrons in the final state [9]), 
they can be determined accurately using present lattice methods. Currently all realistic lat- 
tice calculations of fn^^^ and /^^^j that include the effects of three light dynamical quarks use 
staggered lattice fermions [10, 11] for the up, down, and strange quarks. Because staggered 
fermions are computationally cheaper than other lattice fermion formulations, they allow for 
QCD simulations with dynamical quarks as light as O.OSm^, several lattice spacings, down 
to a ~ 0.045 fm, large physical volumes, and high statistics. This enables lattice deter- 
minations of many light-light and heavy-light meson quantities with controlled systematic 
uncertainties. The results of staggered lattice calculations are largely in excellent numerical 
agreement with experimental results [9]. This includes both postdictions, such as the pion 
and kaon decay constants [12], and predictions, as in the case of the Be meson mass [13]. 
Such successes give confidence that further calculations using the same methods are reliable. 
This is essential if lattice QCD calculations of hadronic weak matrix elements are to be used 
to test the Standard Model and search for new physics. 

The staggered dynamical quark simulations used here employ the fourth-root procedure 
( "rooting" ) for eliminating unwanted extra quark degrees of freedom that arise from lattice 
fermion doubling. The rooting method is not standard quantum field theory, and at nonzero 
lattice spacing it leads to violations of unitarity [14-17] that can be considered nonlocal [18]. 
Nevertheless, there are strong arguments [19, 20] that the desired local, unitary theory of 
QCD is reproduced by the rooted staggered lattice theory in the continuum hmit. Further, 
one can show [15, 21] that the unitarity-violating lattice artifacts in the pseudo-Goldstone 
boson sector can be described, and hence removed, using rooted staggered chiral perturbation 
theory (rSXPT), which is a low-energy effective description of the rooted staggered lattice 
theory [22-24]. When coupled with other analytical and numerical evidence (see Refs. [25- 
28] for reviews and Ref. [29] for a recent study), this gives us confidence that the rooting 
procedure is valid. Indeed, the validity of the rooted staggered lattice simulations is of critical 
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importance to flavor physics phenomenology, since a majority of the unquenched, three-flavor 
lattice results for hadronic weak matrix elements used to determine CKM matrix elements 
and as inputs to constraints on the CKM unitarity triangle come from such simulations [30] . 

In this paper, we present new results for the leptonic decay constants of heavy-light 
mesons containing bottom and charm quarks. We use the "2+1" flavor asqtad- improved 
gauge conflgurations made publicly-available by the MILC Collaboration [31]. These en- 
sembles include the effects of three light, dynamical sea-quark flavors: one with mass rrih 
near (the physical strange-quark mass) and the other two with mass m/ as small as O.lmh. 
We generate light valence quarks for the B and D mesons using the same staggered action as 
in the sea sector, and generate heavy bottom and charm quarks using the clover action [32] 
with the Fermilab interpretation [33]. Because the Fermilab method uses knowledge of the 
heavy-quark limit of QCD to systematically eliminate heavy-quark discretization errors, ex- 
ploiting ideas of Symanzik [34, 35] and of heavy-quark effective theory (HQET) [36-38], it 
is well-suited for both bottom and charm quarks. We simulate with many values for the 
light up/down quark mass (the mass of our lightest pion in both the sea and valence sectors 
is ~ 250 MeV), and at three lattice spacings ranging from a ^ 0.09 fm to a ^ 0.15 fm. 
We then extrapolate our numerical lattice data to the physical up and down quark masses 
and continuum guided by expressions derived in staggered chiral perturbation theory for 
heavy-light mesons (HMSxPT) [39-41]. 

We renormalize the heavy-light axial current with a mostly nonperturbative approach, 
computing the flavor-diagonal (heavy-heavy and light-light) renormalization factors nonper- 
turbatively and then calculate the remaining flavor off-diagonal correction factor [Pa^ ) 
lattice perturbation theory [37, 42, 43]. This procedure has the advantage that is close 
to unity. Furthermore, tadpole diagrams cancel in the ratio needed to obtain , thereby 

Qq 

improving the convergence of the perturbative series. Empirically, the size of the 1-loop 
contribution to Pa4, is found to be small. 

Our results for the charmed-meson decay constants improve upon our published results 
for fi) and fo^ in Ref. [44] in several ways. The coarsest lattices used in this work have a 
smaller lattice spacing (a ~ 0.15 fm) than those used in our previous work (a ^ 0.18 fm). 
The number of conflgurations in the two most chiral ensembles with a 0.12 fm has been 
increased, approximately by factors of 1.4 (sea mi = O.lrrih) and 1.7 (sea mi = O.lAmh)- We 
have added new data on a new a ^ 0.09 fm sea-quark ensemble with a light quark mass 
of O.lmh- We now obtain our results from a combined analysis of our entire data set (all 
partially-quenched mass combinations and lattice spacings). Furthermore, we now compute 
the bottom meson decay constants fs and fs^- We have presented reports on this project 
at several conferences [45-49]; in our flnal analysis of this data set we also improve upon 
bottom and charm quark mass-tuning, with increased statistics and a more sophisticated 
analysis of heavy-quark discretization effects. 

This paper is organized as follows. In Sec. II, we present an overview of the calculation, 
including the gluon and light-quark actions used in generating the gauge conflgurations and 
the light- and heavy-quark actions used in constructing the heavy-light meson correlators. 
We also introduce the mostly nonperturbative method for matching the lattice heavy-light 
current to the continuum, and the treatment of heavy-quark discretization errors from the 
Fermilab action within our chiral-continuum extrapolation. Next, in Sec. Ill, we describe the 
details of our numerical simulations and we present the parameters used, such as the light- 
quark masses and lattice spacings. We also describe the procedure for tuning the hopping 
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TABLE I. Flavor content of the axial vector current and associated CKM matrix element. 



H 
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D 






Ds 




V* 

' cs 


B 




Vub 


Bs 







parameter in the clover action so that it corresponds to b and c quarks. In Sec. IV, we define 
the two-point correlation functions used to extract the decay constant at each value of the 
light-quark mass and lattice spacing. We use two different fitting procedures to obtain the 
decay constants that differ in their treatment of the statistical errors, choice of fit ranges 
and number of states, and choice of input correlators. We include the difference between the 
two in our estimate of the fitting systematic uncertainty. Next, we present the numerical 
details of the calculation of the heavy-light axial-current renormalization factor in Section V. 
Putting the results of the two previous sections together, in Sec. VI, we extrapolate the 
renormalized decay constant data at unphysical quark masses and nonzero lattice spacing 
to the physical light quark masses and zero lattice spacing using HMSxPT. In Sec. VII, we 
estimate the contributions of the various systematic uncertainties to the decay constants, 
discussing each item in our error budget separately. We present the final results for the 
decay constants in Sec. VIII, and compare them to other lattice QCD calculations and to 
experiment. We describe the impact of our results for current fiavor physics phenomenology 
and then conclude by discussing the ongoing improvements to our calculations, and their 
future impact on searches for new physics in the quark fiavor sector. 

Appendix A applies HQET to the Fermilab action to obtain explicit expressions for heavy- 
quark discretization effects. Appendix B contains the complete set of fit results for the heavy- 
light pseudoscalar meson mass and renormalized decay constant for all combinations of sea- 
quark mass, light valence-quark mass, and heavy-quark mass used in the chiral-continuum 
extrapolation. These results will be included as an EPAPS attachment upon publication. 



II. METHODOLOGY 

The decay rate for a charged pseudoscalar meson H (with fiavor content Q and q) to 
leptons is, in the Standard Model, 



.2 



nH^iu) = -^ff,\G^V*^m,ni--^ , (2.1) 



where Mh is the mass of the meson H, Gp is the Fermi constant, and Vqq is the pertinent 
element of the CKM matrix. The decay constant fn parameterizes the pseudoscalar-to- 
vacuum matrix element of the axial vector current, 

{0\A^\H{p))=zp^fH, (2.2) 

where p'^ is the 4-momentum of the pseudoscalar meson. The fiavor contents of the associated 
vector current and CKM matrix element are given in Table I. Note that the neutral Bg 
decays to a charged lepton pair with an amphtude proportional to fs/, hence the CKM 
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factor in the decay rate involves more than one CKM matrix element. Because this process 
is loop-suppressed in the Standard Model, it is potentially sensitive to new physics effects. 
These formulas hold for all pseudoscalar mesons; in the normalization convention used here, 
A(|Kd|/0.97425) = 130.41 ± 0.20 MeV [50]. 

In Eq. (2.2), the 1-particle state assumes the relativistic normalization convention. For 
mesons containing a heavy quark, however, it is more convenient to pull out factors of Mh 
to ensure a smooth Mh — c>o limit: 

(0|^1iJ(p)) {Mh)'^'^ = t{p^/MH)<i>H. (2.3) 

In lattice QCD, the normalization of states on the left-hand side falls out of correlation 
functions more naturally. Thus, most of our analysis, including error analysis, focuses otk^h- 
We then obtain = (pn/V using the experimentally measured value of the meson 
mass [51]. 

To compute the decay constants with lattice gauge theory, we must choose a discretization 
for the heavy quark, the light quark, and the gluons. As in previous work [44, 52-55], we 
choose the Fermilab method for heavy quarks [33] and staggered quarks with the asqtad 
action [56] for the light (valence) quark. The gauge action is Symanzik improved, with 
couplings chosen to remove order a^a^ errors from gluon loops [57], but not those from 
quark loops [58] (which became available only after the gauge-field generation was well 
underway) . 

For heavy bottom and charm quarks, we use the Sheikholeslami-Wohlert (SW) clover 
action [32] with the Fermilab interpretation [33], which connects to the continuum limit as 
aniQ —7- 0. This is an extension of the Wilson action [59], which retains the Wilson action's 
smooth limit as amq — ?■ oo and also remains well behaved for mqa ^ 1. Because this lattice 
action respects heavy-quark spin-flavor symmetry, one can apply HQET to organize the 
discretization effects. In essence, one uses HQET to develop the l/mg expansion both for 
continuum QCD and for lattice gauge theory (LGT) [36-38]. Discretization effects are then 
captured order-by-order in the heavy-quark expansion by the difference of the short- distance 
coefficients in the descriptions of QCD and LGT. Thus, in principle, the lattice heavy-quark 
action can be improved to arbitrarily high orders in l/mq by adjusting a sufficiently large 
number of parameters in the lattice action. (See Ref. [60] for details at dimension 6 and 7. In 
principle, the adjustment can be done nonperturbatively, such as in the scheme of Ref. [61].) 
In practice, we tune the hopping parameter k and the clover coefficient csw of the SW 
action, to remove discretization effects through order l/mq in the heavy-quark expansion. 

The HQET analysis of cutoff effects could be applied to any lattice action with heavy- 
quark symmetry, such as the action of lattice NRQCD [62]. In the latter case, it is simply 
a different perspective on the usual approach to lattice NRQCD, which derives the heavy- 
quark Lagrangian formally, and then replaces derivatives with difference operators. A key 
feature of the Wilson, SW, Fermilab and OK [60] actions is their well-behaved continuum 
limit, which is especially important for charm. For rriqa < 1, one can analyze the cutoff 
effects in a complementary way with the Symanzik effective action [34, 35]. This two-pronged 
attack shows that the difference of short-distance coefficients, mentioned above, vanishes as 
a suitable power of lattice spacing a. In this paper, we shall use our knowledge of this 
behavior to constrain heavy-quark discretization effects in several steps of our analysis. See 
Sees. HI B, VI, and Appendix A for details. 
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The lattice and continuum currents are related by a matching factor Zap^ [37]: 



Za.A^ = + O (a,aA/,(mQa)) + O [o" f ,{mQa)) , (2.4) 

where = denotes equality of matrix elements, and the functions /jj that depend on mqa 
stem from the difference in the HQET short- distance coefficients. In the Fermilab method, 
they remain of order 1 for all values of mqa [33, 60], and they are given explicitly in 
Appendix A. In this work, we compute mostly nonperturbatively [42] and partly in 
one-loop perturbation theory. As shown in the analysis of Ref . [37] , many of the Feynman 
diagrams in the perturbative expansion of Z^^^ are common or similar to those in the flavor- 
conserving renormalization factors ^y*^ and ZyA^ , which can be computed nonperturbatively. 
Therefore, we define p^4 by 



'Z' = P \ Zv4, , (2.5) 

evaluating only p^^ in lattice perturbation theory. 

The flavor-conserving factors account for most of the value of the heavy-light renormal- 
ization factor . They are obtained by enforcing the normalization condition, at zero 
momentum transfer, 

l = Zy.^{H,\V^\\H,), (2.6) 

where Hq is a hadron containing a single quark of flavor g, and V^^ is the lattice version of the 
degenerate vector current. This condition holds for all discretizations and quark masses and, 
hence, the heavy quark {i.e., Zyi^) as well. The remaining correction factor p^^^ is close to 
unity due to the cancellation of most of the radiative corrections including tadpole graphs. 
Although such cancellations have only been explicitly shown at 1-loop in lattice perturbation 
theory [37, 43], we expect similar cancellations to persist at higher orders. Therefore, the 
perturbative truncation error in the heavy-light renormalization factor is subdominant. 



III. LATTICE SIMULATION DETAILS 
A. Parameters 

Table II lists the subset of the ensembles of lattice gauge fields generated by the MILC 
Collaboration [28] used in this analysis. We now describe each entry in the table. 

We analyze data at three lattice spacings: a ~ 0.15 fm, a ^ 0.12 fm, and a ^ 0.09 fm. The 
ensembles contain 2+1 flavors of sea quarks, using the asqtad-improved staggered action [56], 
and the square (fourth) root of the staggered determinant for the two degenerate light sea 
quarks (one strange sea quark). The sea contains one flavor with mass nih close to the 
physical strange quark mass and two degenerate lighter flavors of mass mi. The tadpole 
improvement factor uq is a parameter of the gauge and asqtad staggered (sea) quark action 
and is determined from the fourth root of the average plaquette. We calculate the two-point 
correlation functions on each ensemble from an average over four different time sources. 

The relative lattice scale is determined by calculating ri/a on each ensemble, where ri 
is related to the force between static quarks, r^F(ri) = 1.0 [63, 64]. Table II lists ri/a 
values for each of the ensembles that result from fitting the calculated ri/a to a smooth 
function [65], as explained in Eqs. (115) and (116) of Ref. [28]. 
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TABLE II. The MILC three-flavor lattices and valence asqtad quark masses used in this work. All 
of the valence masses were used in version II of the correlator fits (Sec IV C), while only the ones 
in bold print were used in version I (Sec IV B). 

!v a [fm] aruh ami ri/a riconf X ?^src valence amg 

0.09 0.031 0.0031 0.8779 3.69 435x4 0.0031,0.0037,0.0042,0.0044,0.0052,0.0062, 

0.0087, 0.0124, 0.0186, 0.0272, 0.031 

0.0062 0.8782 3.70 557x4 0.0031,0.0037,0.0044,0.0052,0.0062, 

0.0087, 0.0124, 0.0186, 0.0272, 0.031 

0.0124 0.8788 3.72 518 x 4 0.0031,0.0042,0.0062,0.0087,0.0124, 

0.0186,0.0272,0.031 

0.12 0.05 0.005 0.8678 2.64 678x4 0.005,0.006,0.007,0.0084,0.01,0.012,0.014, 

0.017, 0.02, 0.024, 0.03, 0.0415 
0.007 0.8678 2.63 833x4 0.005,0.006,0.007,0.0084,0.01,0.012,0.014, 

0.017, 0.02, 0.024, 0.03, 0.0415 
0.01 0.8677 2.62 592x4 0.005,0.006,0.007,0.0084,0.01,0.012,0.014, 

0.017, 0.02, 0.024, 0.03, 0.0415 
0.02 0.8688 2.65 460x4 0.005,0.006,0.007,0.0084,0.01,0.012,0.014, 

0.017, 0.02, 0.024, 0.03, 0.0415 
0.03 0.8696 2.66 549x4 0.005,0.006,0.007,0.0084,0.01,0.012,0.014, 

0.017, 0.02, 0.024, 0.03, 0.0415 
0.15 0.0484 0.0097 0.8604 2.13 631x4 0.0048,0.007,0.0097,0.013,0.0194,0.0242, 

0.029,0.0387,0.0484 
0.0194 0.8609 2.13 631x4 0.0048,0.007,0.0097,0.013,0.0194,0.0242, 

0.029,0.0387,0.0484 
0.029 0.8614 2.13 576x4 0.0048,0.007,0.0097,0.013,0.0194,0.0242, 

0.029,0.0387,0.0484 



In order to fix the absolute lattice scale, one must compute a physical quantity which can 
be compared directly to experiment. The combination of the PDG's value of /,r with MILC's 
2009 determination of rif^^ [66] yields ri = 0.3117(6) (^3^) fm. From an average of three 
methods for scale setting, including one based on T splittings, the HPQCD collaboration 
obtains ri = 0.3133(23) (3) fm [67], consistent with MILC. Symmetrizing MILC's error range 
gives ri = 0.3108(21) fm, and a straightforward average with the HPQCD result then yields 
ri = 0.3120(16) fm. This average omits likely correlations, due to the use of MILC sea-quark 
configurations by both groups. Conservatively assuming a 100% correlation, we inflate the 
error to 0.0022 fm. Finally, for convenience, we also shift the central value slightly, back to 
the 2009 MILC central value. We thus take n = 0.3117(22) fm in this paper. 

The complete list of light (asqtad) valence quark masses simulated in this analysis is 
also given in Table II. The mass values are selected to be roughly logarithmically spaced, but 
to also include the set of light sea quark masses simulated at each lattice spacing. We use a 
multimass solver to compute the valence quark propagators. The marginal numerical cost 
of including masses heavier than our lightest rriq ~ O.lm^ is small and logarithmic spacing 
is designed to constrain the chiral logarithms. 
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TABLE III. Table of clover coefficients and n values for charm and bottom used in heavy-light 
two-point simulations. 



~ a [fm] 


ami/amfi 




charm bottom 


0.09 


0.0031/0.031 


1.478 


0.127 0.0923 




0.0062/0.031 


1.476 






0.0124/0.031 


1.473 




0.12 


0.005/0.05 


1.72 


0.122 0.086 




0.007/0.05 


1.72 






0.01/0.05 


1.72 






0.02/0.05 


1.72 






0.03/0.05 


1.72 




0.15 


0.0097/0.0484 


1.570 


0.122 0.076 




0.0194/0.0484 


1.567 






0.0290/0.0484 


1.565 





In Table III, we show the coefficient of the Sheikholeslami-Wohlert term csw of the clover 
action and the k values used to compute heavy-light two-point functions. The coefficient of 
the clover term is set to the tadpole-improved tree-level value csw = Uq^- For the a ~ 0.09 
and 0.15 ensembles the tadpole coefficient is taken from the average plaquette. We note, 
however, that at lattice spacing a ~ 0.12 fm the tadpole coefficient Uq appearing in both the 
valence asqtad action and the heavy quark clover action is taken from the average of the 
Landau link evaluated on the ami/am^ = 0.01/0.05 ensemble. Hence, in our a ~ 0.12 fm 
lattice data there is a mismatch between light valence and sea quark mass definitions. As 
discussed in Sec. VII, this (inadvertent) choice leads to a small error in the decay constants. 
We have remedied this mismatch by using the plaquette uo everywhere in new runs started 
while this analysis was underway. 

The charm and bottom kappa values listed in Table III are based on our initial kappa 
tuning analysis using about one fourth of our final statistics. We then used a larger data set 
to refine our determination of the k values corresponding to bottom and charm as described 
in the next subsection. We adjust our data post-facto to correspond to tuned values of Kc 
and Kb using the measured value of the derivative 6(f)/ 6k. 

B. Input quark masses rric and mf, 

Our method for tuning k, for charm and bottom quarks closely follows that of Ref . [55] , 
where further details can be found. We start with the dispersion relation for a heavy particle 
on the lattice [33] 

E'ip) = Ml + + \a, {ap'f + lA,a' 1^/ + • • • ' (3-1) 

where 

Ml = E{0) (3.2) 
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TABLE IV. Hopping-parameter values used to compute the dispersion relation. 



~ a [fm] 


^conf ^ IT-src 


charm 


bottom 




0.09 


1912 X 4 


0.1240, 0.1255, 0.1270 


0.090, 0.092, 


0.094 


0.12 


592 X 4 


0.114, 0.117, 0.119, 0.122, 0.124 


0.074, 0.086, 


0.093, 0.106 


0.15 


631 X 8 


0.100, 0.115, 0.122, 0.125 


0.070, 0.076, 


0.080, 0.090 



is called the rest mass, and the kinetic mass is given by 

-1 _ dE{p) 



(3.3) 

p=0 



These meson masses differ from corresponding quark masses, mi and m2, by binding-energy 
effects. The bare mass or, equivalently, the hopping parameter k must be adjusted so that 
these masses reproduce an experimental charmed or 6-flavored meson mass. When they 
differ, as they do when mga ^ 1, one must choose. Decay constants are unaffected by 
the heavy-quark rest mass mi [36], so it does not make sense to adjust the bare mass to 
Ml. We therefore focus on M2, adjusting k to the strange pseudoscalars Dg and B^, both 
because the signal degrades for lighter spectator masses and because this avoids introducing 
an unnecessary systematic uncertainty due to a chiral extrapolation. 

The first step is to compute the correlator C2 ^ ^ {t,p) in Eq. (4.8) (below) for several 
3-momenta p and several values of k, and light quark mass, bracketing charm and bottom, 
and strange, respectively. We use all momenta such that |p| < Att/L. Second, we fit the 
time dependence of the multichannel correlation matrix ^ ^ to a sum of exponentials — 
including the usual staggered-fermion oscillating terms — and extract the ground state energy 
aE{p) by minimization of an augmented [55, 68, 69]. Third, we fit the energies to the 
dispersion relation given in Eq. (3.1), through 0{pf). The output of this fit is aMi, M1/M2, 
A4, and 744/, all as functions of k. Fourth, we form M2(/«) from the first two fit outputs 
and Ti/a, propagating the error with a single-elimination jackknife. Finally, we obtain our 
tuned Kc and Kb by interpolating in k so that M2{k) matches the experimentally known 
Ds and Bg masses. The n values used to compute M2 are listed in Table IV. For each of 
the lattice spacings listed, we used the ensemble with light-to-strange sea-quark mass ratio 
ami/anih = 0.2. The resulting tuned values of Hc and Kb are shown with errors in Table V. 

We constrain the coefficients A4 and A4/ with Gaussian priors derived from the HQET 
theory of cutoff effects, adding the contribution of the priors to the in the minimization 
procedure [68, 69]. (In principle, we could include such priors for Mi and M1/M2 too, but in 
practice we take priors so wide that these fit parameters are solely data-driven.) Neglecting 
binding energies, we have exact tree- level expressions for 04 and 04/, the quark analogs 
of Ai and A^i. The differences A4 — 04°^ and A4/ — a!^} stem from both perturbative and 
nonperturbative effects. The asymptotics of the former can be estimated along the lines of 
Appendix A 3, and the latter can be deduced following the methods of Refs. [36, 70]. Briefly, 
we constrain y4„(/t), n G {4,4'}, to a Gaussian with central value 

a|^'(moa) + asa\^\moa) + AaA'^{moa). (3.4) 

Here is the exact tree-level contribution, a™' is an estimate of the one-loop contribution, 
and A'^ is an expression for the binding-energy contribution. The width of the Gaussian is 
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TABLE V. Hopping parameter values Kc and Kb corresponding to charm and bottom. The outputs 
of the tuning are labeled Ktunedj where the first error is from statistics and the second is from ri, 
which enters through matching to the experimentally-measured Dg and Bg meson masses. The 
derivative dcp/dn is used to correct the values of (p obtained with the simulated values Kgim listed 
in Table III to the tuned values given below. 





charm 




bottom 




~ a [fm] 




d(f)/dK 




d(t)/dK 


0.09 


0.12691(18)(13) 


-21.66 


0.0959(13)(3) 


-7.41 


0.12 


0.12136(37)(19) 


-18.23 


0.0856(19)(3) 


-6.82 


0.15 


0.12093(36)(24) 


-15.40 


0.0788(11)(3) 


-6.07 



determined by combining in quadrature the chosen widths of the separate contributions, as 
outhned in Appendix A3. 

The details of the and Kb determination differ from that of Ref. [55] in two respects. 
First, we use the pseudoscalar meson masses rather than the spin average of pseudoscalar 
and vector meson masses, leading to a modest reduction of the statistical error. Second, 
we include the quartic terms in Eq. (3.1), allowing us to fold discretization effects directly 
into the dispersion-relation fit. Although we consider these two changes improvements, the 
change in the tuned k values as compared to Ref. [55] stems primarily from the substantial 
increase in statistics on key ensembles. 

The computations of the correlation functions needed to extract (po and 0b have been 
carried out using the fiducial values listed in Table III. These simulation k's were obtained 
near the beginning of the project, but while the runs were in progress, we redetermined the 
hopping parameters utilizing increased statistics and reflecting an updated value of ri [66]. 
The resulting improved determinations of Kc and Kb differ slightly from the simulation values. 
In order to adjust from the simulated value Kgim to the tuned value Ktuned, we write 

tuned 0sim ~l~ ('^tuned ^sim)? ('^■^) 

where the derivatives d(f)/dn listed in Table V are obtained from tuning runs with nearby k 
values. As explained in Sec. VII, these derivatives are also used to propagate to the decay 
constants the statistical and scale uncertainties on Ktuned listed in Table V. 



IV. TWO-POINT CORRELATOR FITS 

We obtain the unrenormalized decay amplitude for every combination of heavy-quark 
mass, light-quark mass, and sea-quark ensemble by fitting the heavy-light meson two-point 
correlation functions, described in Sec. IV A. We use two independent fitting procedures, 
which we refer to as "Analysis I" and "Analysis II". These procedures differ in several re- 
spects. In Analysis I, we use a jackknife procedure for estimating errors, while in Analysis II, 
we use a bootstrap procedure. The two analyses also differ in their methods for handling 
autocorrelations in the data and in their choices of fit ranges, priors for masses and ampli- 
tudes, and numbers of states included. In the end, we use Analysis I (Sec. IV B) to obtain 
central values, and use differences from fits with different distance ranges and/or number 
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of states included, and from Analysis II (Sec. IV C) to estimate the systematic error due to 
choices made in the fit procedure. 



A. Lattice correlators 



The lattice axial-vector current is given by 

Atix) = [^{xh'^'n{x)U{x), (4.1) 

where xi^) is the one-component field appearing in the staggered action, and Q{x) = 
^xi/a^x2/a^xs/a^x4/a transformation connecting naive and staggered fields [71]. The 

heavy-quark field \1/ is a four-component (Dirac) spinor field, and the remaining free Dirac 
index is interpreted as a taste label. 

To remove tree-level discretization errors in the lattice axial current, the heavy-quark 
field ^ is "rotated" : 

^ = [1 + adi7 ■ D]^, (4.2) 
where is the field appearing in the clover action. Tree-level improvement is obtained when 

di = 7 T, 4.3 

2 + moa 2(l + moa) ^ ^ 

where 

moa= — (—-- ) (4.4) 

is the tapdole-improved bare mass. The critical hopping parameter Kent is the one for which 
the clover-clover pion mass vanishes. 

As usual in lattice gauge theory, we obtain the matrix element in (2.3) from two-point 
correlation functions. We introduce pseudoscalar operators 

Oi'\x) = J2[^{y)S{y,xh'n{x)lxixl (4.5) 
y 

depending on a "smearing" function S. In this work, we use two functions, the local (or 
unsmeared) source S{x,y) = Sxy, and the smeared source (in Coulomb gauge) 

S{x,y) = 6a:^y^S{x - y), (4.6) 

where S{r) is the IS solution of the Richardson potential for the quarkonium systems [72]. 
We obtain S{x — y) by scaling the radial Richardson wavefunction to lattice units, in- 
terpolating it to lattice sites, and then using it as the spatial source for the heavy-quark 
propagators [73]. 

We introduce two-point correlation functions 



^fW = EEK(^'^)^i'nO)), (4.7) 

a=l X 
4 

Ci'^'^\t,p) = Y,J2^'''-'' (Oi'^^\t,x)0i'^\0)) , (4.8) 



a=l X 
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where (•) now represents the ensemble average. For large time separations, $2 is propor- 
tional to the matrix element 0^, and the proportionality is determined from C2 it,0). 
Each two-point function is constructed from a staggered quark propagator with local (S) 
sources and sinks. We compute C2 functions for all (four) combinations 5*1, 5*2 = S and 
IS, requiring heavy clover quark propagators with all combinations of IS smeared and local 
sources and sinks. Only the local sink clover propagators are needed to compute the $2 
functions. With the sum over tastes in Eqs. (4.7) and (4.8), the correlation functions $2 
and C2 can also be cast in a heavy-naive formalism [74]. 

The staggered light quarks in the axial-current and pseudoscalar two-point correlation 
functions lead to the presence of opposite-parity states that oscillate in time as (—1)*- Hence 
the two-point functions take the following forms: 

$f )(t) = ^Sf) (e-*^* + e-*^(^-*)) + Ai'^ {-if (e-*^* + e'^^^^"*)^ 

+ Af)(e-*^'* + e-*^'(^-*))+..., (4.9) 
^^^)(t,p = 0) = (e-^^* + e-*^(^-*)) + 1(^1)1(^2) (-1)* (^e-^* + e-*^(^-*)) 

where a prime denotes a standard excited state of the same parity and a tilde denotes the 
mass or amplitude of an opposite-parity state. The oscillating behavior is visible throughout 
the entire lattice temporal extent, and must be included in fits to extract the ground-state 
mass and amplitudes. 

We then obtain the renormalized decay amplitude in lattice units from the ratio 

a'^'^H = ^^J(s^^ (4-11) 

where and A^^'^ are the amplitudes of the ground state exponentials defined in Eqs. (4.9) 
and (4.10), and the renormalization factor is discussed in Sec. V. 



B. Analysis I 

Our primary analysis of two-point correlation functions $2 C2 ^ ^ — "Analysis I" — 

proceeds as follows. The amplitudes and A'^^'' in Eq. (4.11) are determined from fits 
to multiple correlators using the full data correlation matrix. In Analysis I, we start by 
fitting combinations A, B, C and D in Table VI. We find combination A, which uses the 
axial-current correlator with a IS" smeared source and the pseudoscalar correlator with a 
IS" smeared source and sink, to be suitable. The extra complexity of combinations of three 
correlators (C and D) give little benefit, and the errors from combination A are somewhat 
smaller than those from combination B. 

For fits to charm-light meson correlators, we include just one simple exponential (the 
desired state) and one oscillating exponential, which we call a "1-1-1 state fit". We choose 
the minimum distance, tmm, such that contributions from excited states are small compared 
to our statistical errors. Because we fit two propagators simultaneously while imposing the 
constraint that the masses be equal, this is a six parameter fit: two amplitudes for each 
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TABLE VI. Combinations of two-point functions that can be used to extract a?/'^(f)H- All combina- 
tions of two and three correlators are shown. Additional combinations of four or more correlators 
are not enumerated. 



two-point 


fit combinat 


ion 


function 


A B C D E 


F 




• • • • 


• 










• • 


• 




• • 






• 


• 




• 





propagator and a common mass for each of the simple and oscillating exponentials. To help 
stabilize the fit, the amplitudes and mass of the oscillating state are weakly constrained by 
Gaussian priors, which are incorporated as additional terms in in the fitting procedure [68, 
69] . The central values for these priors are determined by a trial fit where the prior for the 
opposite parity mass is set to 500 ±250 MeV above the ground state mass and the amplitudes 
are unconstrained. Then the jackknife fits use central values for the opposite parity state 
amplitudes and mass determined by the trial fit, with widths (Gaussian) that are typically 
three to ten times the error estimates on these parameters, so in the end the priors make 
a negligible contribution to x^- (Although 500 MeV is a reasonable guess for the mass gap 
to the first excited state of the meson, we actually expect that this excited state in the fit 
approximates the contributions of a number of physical states, likely including both single 
and multiparticle channels.) Empirically, the width of the prior is made narrow enough 
to insure that the fits converge to sensible values. We propagate the uncertainties in the 
correlator fits to the subsequent chiral- continuum extrapolation with a jackknife procedure. 
In the jackknife resamples, we center the priors at the values found in the fit to the full 
ensemble, again with widths that are typically three to ten times the error estimates on 
these parameters. 

The bottom meson correlators fall off much more rapidly with t, so it is difficult to take 
a large enough minimum distance to insure that excited state contributions are negligible. 
Therefore we use a fit with two simple exponentials and one oscillating exponential or a 
"2+1 state fit". The mass of this excited state is also weakly constrained by priors in the 
same way that the opposite parity mass is, except that the width of the prior on the excited 
state mass is set to 200 MeV. 

Figure 1 shows the heavy-light pseudoscalar mass as a function of the minimum time used 
in the fit. The left-hand plots show sample fits to bottom correlators, while the right-hand 
plots show sample fits to charm correlators. We select fitting ranges to give reasonable fits 
for all sea-quark ensembles and all valence-quark masses. We quantify the goodness-of-fit 
with the "p value" [51], which is the probability that a fit with this number of degrees of 
freedom would have a larger than this value. Table VII gives the fit ranges for charm- 
light and bottom-light correlators on the three lattice spacings, both for the fits used for 
the central values and for alternate fits used in estimating systematic errors from choices 
of fit parameters. The meson masses, a^^'^cpH values, and p values for the data set used in 
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TABLE VII. Numbers of states and time ranges used in two-point Analysis I. In the number of 
states, "1+1" means one simple exponential and one oscillating state (opposite parity). The fits 
in columns two through five were used for the central values, while the fits in columns six through 
nine were used in estimating systematic errors from the choice of fit ranges (see Sec. VII C). 





central fits 


alternate fits 




charm 


bottom 


charm 


bottom 


~ a [fm] 


^states t range 


?T-states t range 


^states t range 


nstates t range 


0.15 


1+1 11-23 


2+1 4-20 


1+1 12-23 


1+1 9-20 


0.12 


1+1 14-31 


2+1 5-22 


1+1 16-31 


1+1 12-22 


0.09 


1+1 21-47 


2+1 7-30 


1+1 23-47 


1+1 16-30 



Analysis I are tabulated in Appendix B. 

The decay amplitude a^^'^cpn is highly correlated among different light valence-quark 
masses on the same ensemble. To propagate the correlations among the different va- 
lence masses to the subsequent chiral-continuum extrapolation, in Sec. VI, we use a single- 
elimination jackknife procedure to estimate the covariance matrix of a^^'^cpn for the selected 
valence quark masses. This is done by computing the covariance matrix of the single elimina- 
tion jackknife samples, and multiplying by (A^ — 1)^, where N is the number of configurations 
in the ensemble. In fact, when all valence quark masses are kept, the covariance matrices 
are close enough to singular to be unmanageable. This reflects the fact that the correlators 
for intermediate valence masses can be very accurately predicted from the correlators for 
nearby masses, so some of the correlators provide very little new information. Therefore, we 
omit some valence quark masses, using only those set in bold in Table II. 

We use a single elimination jackknife rather than an omit-J jackknife because a large 
number of samples are needed to compute a reliable covariance matrix. Successive config- 
urations in the ensemble are not independent, however, so we must take autocorrelations 
into account. We do so by repeating the calculation after first blocking the data by a factor 
of four. (This block size of four is determined from tests on the a ~ 0.12 fm lattices using 
fit Analysis II, for which it gives a reasonable compromise between suppressing autocorre- 
lations and leaving enough data points for the statistical analysis.) We then compute, for 
each valence-quark mass i, the ratio Ri of the diagonal element of the covariance matrix 
with a block size of four to the same element of the unblocked covariance matrix: 

= alP/au, (4.12) 

where the superscript denotes the block size. The rescaled covariance matrix for a?/'^(j)H is 
given by 

C\f = C^j^R^R,, (4.13) 

which preserves the eigenvalue structure of the covariance matrix, whereas simply using 
the covariance matrix of the blocked data would be more likely to produce spurious small 
eigenvalues. The rescaling factors Ri themselves have errors, and in many cases turn out 
to be less than one. In such cases, we do not replace the Ri by one, despite the fact 
that this would likely be a better estimate of the individual Ri. Doing so would yield a 
covariance matrix with a bias toward larger errors, and could produce misleading estimates 
of goodness-of-fit in the later analysis. 
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FIG. 1. Ground-state rest mass Mh versus minimum distance tmin included in the fit. For each 
lattice spacing, we show an ensemble with the dynamical light-quark mass mi in the middle of the 
range. Similarly, we show correlators with a valence quark mass my in the center of the ranges 
used. The top two panels are at a ~ 0.15 fm, the middle two at a ~ 0.12 fm and the bottom two at 
a ~ 0.09 fm. In each row the left panel shows results for charm and the right-panel shows results 
for bottom. The size of each plot symbol is proportional to the p value (confidence level) of the fit, 
with the symbol size in the legends of the upper right panel corresponding to p = 50%. The red 
octagons are for fits including one state of each parity ("1+1 fits") and the blue squares are for fits 
including an excited state of the same parity as the ground state ("2+1 fits"). In each panel, the 
arrow indicates the fit that is used in Sec. VI. 
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Finally, we combine the covariance matrices from all of the individual ensembles into 
larger covariance matrices, one each for the charm and bottom Since different 

ensembles are statistically independent, these large covariance matrices are block diagonal, 
with each block containing the correlations among different light valence-quark masses on a 
single sea-quark ensemble. 

C. Analysis II 

Analysis II is a second, independent analysis of the two-point correlators that uses the 
bootstrap method to propagate correlated errors from the two-point analysis through to the 
chiral extrapolations. In Analysis II, we block average the two-point correlator data over four 
sequential configurations (which themselves are spaced by more than four trajectories) before 
the analysis. In the bootstrap procedure, we resample the data (with replacement), taking 
the number of sampled configurations to be equal to the number of blocked configurations 
in each bootstrap ensemble. For each bootstrap ensemble, we recompute the covariance 
matrix. During the bootstrap process, we randomly draw from a gaussian distribution new 
prior mean values of each constrained parameter belonging to an excited state while keeping 
the widths fixed. The ground state parameters are given loose priors so that the fitted values 
are determined by the data. To help stabilize the fits, the ground state prior means are not 
randomized in the bootstrap. Prior values for the energy splittings are taken from a chiral 
quark model calculation for the D and B meson systems [75]. Prior widths are taken to 
be about 200 MeV for excited states. Excited state amplitudes log(y4^'^)) have a prior width 

<^\ogA = 2. 

On each gauge ensemble, the same sequence of gauge resamplings is taken for all valence 
niq to preserve correlations among a^/'^cf)^ values. Our final results are based upon 4,000 
bootstrap replications of the data. We use the central values of from the fits to 

the entire ensemble in the chiral-continuum extrapolation, and use the bootstrap values to 
obtain the covariance matrix. 

To optimize the determination of a^^'^cpn, we compare simultaneous fits of up to six two- 
point functions; the various combinations of up to four functions are listed in Table VI. At 
a minimum, one axial-current correlator must be paired with one propagator (combinations 
A or B in Table VI) to extract a?/'^(j)H- Combination A, using smeared operators, is used 
in Analysis I, described above. Because fits of four or more two-point functions over a 
wide time range can lead to a poorly determined data covariance matrix having large rank 
relative to the number of available configurations, we focus on combinations having two or 
three correlators. Unlike combination A, combination B does not take advantage of smeared 
sources and the ratio does not show convincing plateaus over the range of times with decent 
signal to noise. Comparing combination C to D, the smeared source in C is less noisy than 
the smeared sink in D. 

Given these considerations, for fits to charm correlators, we use two-point function com- 
bination C to obtain a^^'^cpn which uses both of the axial current functions. We look for 
stability of the ground-state mass and amplitude when varying tmin, tmax, and the number 
of excited states included in the fit. We also compare fit results from other combinations 
of correlators to check that we have isolated the correct ground-state energy and matrix 
element. Our final results come from fits accounting for two pseudoscalar states and two 
(oscillating) opposite-parity states. 

For fits to bottom correlators, we use combination B for our final results; this is the same 
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set used in Analysis I. Combination C gives fits with rather low confidence levels for the 
B meson and tends to result in larger errors for a^^'^cp^. Again, we examine fits varying 
the time range; we also try fits with up to three pseudoscalar states plus three oscillating 
opposite parity states. We use these fits and fits to alternate combinations of two-point 
correlators consistency check. 

The fit results from the two different analyses are consistent with each other for most 
cases, but there are a few cases where they differ by a standard deviation or more (see 
Figure 11). The a^^'^cpH results from the two analyses are propagated through the chiral- 
continuum extrapolations in Sees. VI B and VI C. The resulting differences in the extrap- 
olated results in turn provide the basis for our systematic error analysis due to fit choices 
given in Sec. Vll C. 



V. HEAVY-LIGHT CURRENT MATCHING 

In this section, we discuss in more detail the ingredients of Eq. (2.5), which allow a 
"mostly nonperturbative" matching procedure [42]. 



A. Perturbative calculation of Pa4 

The perturbative expansion of pAi can be written as 

Qq 



pAi = 1 + as(g*)p||l(mQa,m(,a) + . . . . (5.1) 



where is the strong coupling and pj^l is the one-loop coefficient. The one-loop coefficient 

is calculated in Ref. [43] using lattice perturbation theory, where we see explicitly that p|^i is 
small because most of the one- loop corrections cancel. The self-energy contributions cancel 
exactly (to all orders, in fact), and, in practice, we are in a region where p|^i(mQa, mgo), 
viewed as a function of mga, has two zeroes. Therefore the renormalization factor is 
close to unity for both bottom and charm. 

The perturbative calculation of p^4 in Eq. (5.1) proceeds as follows. We use ^^(g*) 
defined in the V scheme [76] as determined in Ref. [77], and take q* = 2/a, which is close 
to the optimal choice of Refs. [76, 78] for a wide range of quark masses. The one-loop 
coefficients p|^i are computed for light-quark masses arriq = 0.001,0.01,0.04 to cover the 
range used in this analysis. From these we obtain p 44 at other light-quark masses by linear 

Qq 

interpolation in anig. For illustration. Table Vlll lists p^4 and pA^^ evaluated at the light 
valence mass arriq = 0.01 for the eleven sea-quark ensembles used in this work. Note that 
the sea-quark mass dependence is indirect, via the plaquette used to determine as{q*). The 
dependence on the light-quark mass in the current is very mild: for bottom, p^4 changes 

bq 

with aruq by 0.07-0.2%, depending on lattice spacing, and for charm, p^4^ changes by around 
0.1%. On the fine ensembles, the arriq dependence is almost as large as the total one-loop 
correction because the overall cancellation, especially in p^4 , is so fortuitously good. 
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TABLE VIII. The perturbative correction factor /j^^ for the heavy-hght current A'^ at the sim- 
ulated charm and bottom heavy quark k values given in Table III and at amq = 0.01 for the 
different sea-quark ensembles. The statistical errors associated with the numerical integration are 
negligible. 



~ a [fm] 


ami /airih 


bq 




0.09 


0.0031/0.031 


1.0026 


1.0000 




0.0062/0.031 


1.0026 


1.0000 




0.0124/0.031 


1.0026 


1.0000 


0.12 


0.005/0.05 


1.0081 


0.9959 




0.007/0.05 


1.0081 


0.9959 




0.010/0.05 


1.0081 


0.9959 




0.020/0.05 


1.0080 


0.9960 




0.030/0.05 


1.0079 


0.9961 


0.15 


0.0097/0.0484 


1.0270 


0.9937 




0.0194/0.0484 


1.0267 


0.9938 




0.0290/0.0484 


1.0265 


0.9938 



B. Nonperturbative computation of .Z^^4 and ^v^q 

The nonperturbative part of the matching factor Zaa is obtained from the temporal 

Qq 

components of the clover-clover and staggered-staggered vector currents. At zero-momentum 
transfer, the (correctly normalized) vector current simply counts flavor-number, so it is 
possible to obtain Zyi nonperturbatively for any discretization and any mass [42]. 
For the staggered-staggered current, we compute 

Cf^''\h, 0, ti) = E E y)V»0f^^\h, X)) , (5.2) 

ah x,y 

where, as in Eq. (4.5), 0a is a smeared or local clover- staggered meson operator with mass 
chosen to optimize the signal, and 

KK^) = xix)p{xh'n{x)Ux{x) (5.3) 



is the temporal component of the staggered-staggered vector current. The three-point func- 
tions Cs are computed from the same staggered quarks used for the clover- staggered two 
point functions. The staggered quark is transformed into an improved naive quark by apply- 
ing the Q matrix; this naive quark at time ti is then used as the source term when computing 
the charm propagator. We smear the source at ti so that Si = S2. 

We compute Zy4 using a Dg meson [cf. Eq. (4.5)], which provides a good signal. The 

three-point function Cg ^ ^ (t2,0,ti) contains states of both the desired and the opposite 
parity, with the latter carrying oscillating (—1)* dependence. We construct C3 ^ (t2,0,ti) 
with local sources Si = S2 = 6 and compute it at multiple even and odd values of ti and t2 
in order to disentangle the ground-state amplitude from the other contributions. 
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TABLE IX. Light-lieht vector current renormalization factor Zv^ ■ Values in bold are used in 
computing the heavy-light current renormalizations. With our conventions, the tree-level value of 
Zy4^ is 2. A colon is used to represent the range of time values included in the fit. 

a [fm] ami/amh riconi —ti t2 amg Zyi^ 

0.09 0.0124/0.031 518 23:12 11:13 0.0272 1.868(49) 

23:12 11:13 0.0124 1.883(69) 
0.12 0.01/0.05 592 15:9 7:11 0.03 1.853(45) 

0.007/0.05 523 20:7 7:12 0.03 1.882(56) 

0.15 0.0097/0.0484 631 20:5 4:12 0.0484 1.704(34) 

20:5 4:12 0.029 1.709(40) 
20:5 4:12 0.0242 1.711(42) 
20:5 4:12 0.0194 1.707(45) 
20:5 4:12 0.0097 1.662(55) 



Within the time range ti < < t2 and in the hmit of large separations, |ti|,t2 ^ o, 
Ci'''\t2,0,h) = Z-lA'exp{-E{h-ti)) 

99 

+ Z'AB [(-l)*i exp {E'ti - Et2) + (-1)*^ exp {Et^ - E%)\ 

+ Z"52(-l)(*i+*2) exp {-E\t2 - ti)) + . . . , (5.4) 

neglecting contributions from excited states. We extract Zyi^ from a minimum fit to the 
three-point function using the right-hand side of Eq. (5.4) as the model function. The fit is 
linear in the free parameters Zzl , Z" and Z', while we fix the ground-state energies E and 

99 

E\ and the operator overlaps A and B to the values determined by fitting the two-point 
function Cf^(t,0). We use a single-elimination jackknife procedure to compute the data 
covariance matrix. 

Table IX presents our results for ZyA on the ensembles used in this work. The three-point 
functions for the ZyA calculation are generated at a single source time, tsrc = (instead of 

the four used for two-point functions $2 ^2 ^ ^ )■ At a ~ 0.12 fm we have results 

at two values of the sea quark masses which are consistent within errors. At a 0.09 and 
0.15 fm we have results for several values of mq. We do not see evidence for a dependence 
upon TUq with current statistics. The errors, however, increase at smaller quark mass. Hence, 
we use the ZyA^ corresponding to ~ in Eq. (4.11). In the table, they are set in bold. 
For the clover-clover current, we compute 

^^)(t2, ti, 0) = 5^ {&'^^Kh. y)v^Q{h, x)d^'^\o)) , (5.5) 

where 

V^q{x) = ^{x)^'^{x) (5.6) 

is the temporal component of the (rotated) clover-clover vector current. The clover-clover 
bilinear ^ 

0^'\x) = Y,Hy)S{x,y)j's{x) (5.7) 
y 
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TABLE X. Heavy- heavy vector current renormalization factor Zy4^ computed at several k values, 
covering the charm and bottom quark masses, for three lattice spacings. 



a [fm] ami/am^ 



^conf ^ '^src 



0.09 0.0062/0.031 1912 x 2 



0.1283 0.2749(4 
0.127 0.2830(4 
0.110 0.3856(6 
0.0950 0.4730(8 
0.0931 0.4840(9 



0.12 0.007/0.05 2110 x2 



0.124 
0.122 
0.116 
0.098 
0.086 
0.074 



0.2899(4 
0.3028(4 
0.3410(5 
0.4507(7 
0.5209(10) 
0.5894(15) 



0.15 0.0194/0.0484 631 x 2 



0.122 
0.118 
0.088 
0.076 



0.3195(14) 
0.3440(16) 
0.5195(48) 
0.5898(81) 



consists of a heavy-quark field corresponding to charm or bottom, as the case may be, 
and a hght clover-quark field s with mass chosen to provide a good signal. At large time 
separations, these three-point functions are proportional to Z~} , with the proportionality 

coming from 

Ci^'^'\t) = (^d^^'^\t,x)d^'^^\0)^ . (5.8) 

We compute Zy4^ using a Qs meson, where the strange quark is simulated with the 
clover action to circumvent oscillating opposite-parity states [cf. Eq. (5.7)]. We restrict 
our calculation of C2,3 to 5* = 5*1 = 5*2 using both local and IS smearing functions. The 
function C2 combines a local-local clover quark with mass around and a heavy clover 
quark propagator with source and sink S. The function C3 requires the same heavy- and 
light-quark propagators as needed in C2. An additional heavy-quark propagator originating 
from t2 has as its source the hght quark propagator restricted to t2, multiplied by 75 and 
convolved with smearing function S. 

In Eq. (5.7), we use a random color wall source with three dilutions for both the heavy and 
light spectator quarks that originate from t = 0. We generate two- and three-point functions 
for both local-local and smeared- smeared source-sink combinations where the smearing is 
applied to the heavy quark. We compute the 2- and 3-point functions at several values of k 
spanning a range from around the charm quark to the bottom quark. We determine Zy4^ 
from a fit to the plateaus in the jackknifed ratio of the three-point and two-point functions. 
Our results are summarized in Table X. 

In order to properly normalize the derivative d(j)/dnQ (see Eq. (3.5)), we need values 
of Zt/4 at K values other than those used in the Zv* simulations. We therefore fit the 
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TABLE XI. Heavy-heavy vector current renormalization factor '^v^^ corresponding to the charm 
and bottom Ksim values used in the decay constant simulations. 





charm 


bottom 


~ a [fm] 


^QQ 




0.09 


0.127 0.2829(4) 


0.0923 0.4891(9) 


0.12 


0.122 0.3029(4) 


0.086 0.5216(10) 


0.15 


0.122 0.3199(14) 


0.076 0.5868(81) 



o 

CD 

I 



O 

o > 
N 



1.00- 
0.95- 
0.90- 
0.85- 
0.80- 
0.75- 
1.00- 
0.95- 
0.90- 
0.85- 
0.80- 
0.75- 
1.00- 
0.95- 
0.90- 
0.85- 
0.80- 
0.75- 



0.2 



0.4 0.6 
mo/(1 +mo) 



0.8 



1.0 



FIG. 2. Plot of Zyi /(I — Quqk) vs. moa/(l + mga) for the three lattice spacings. 



simulation results to the interpolating quartic polynomial 

4 

Zy.^{K) = l + Y,c,^^' (5.9) 

which reproduces the infinite mass limit Zyi — 1. Our codes employ the hopping parameter 
version of the action; so, at tree level ci = —Quq and for j > 1, Cj = 0. We constrain the 
interpolation parameters to the tree-level values taking aj = 4 as the widths. Table XI 
shows values for Zy4^ interpolated to the nominal charm and bottom Ksim used in our 
decay constant runs. Figure 2 plots the data in Table X together with the interpolation of 
Eq. (5.9). To aid perturbative intuition, the values of Zy4^ in the figure are scaled by the 

tree-level expression 1 — Quqk; the relation between k and moa/(l + rriQa) can be inferred 
from Eq. (4.4). 
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VI. CHIRAL AND CONTINUUM EXTRAPOLATION 



In this section, we present the combined chiral and continuum extrapolations used to 
obtain the physical values of the and D(^s) meson decay constants. We first discuss 
the use of SU{3) chiral perturbation theory for heavy-light mesons in Sec. VIA, giving the 
formulas used for the chiral fits and describing our method for incorporating heavy-quark 
and light-quark discretization effects into the extrapolation. We then show the chiral fits for 
the D system in Sec. VI B, and for the B system in Sec. VI C. 



A. Chiral Perturbation Theory framework 

The errors introduced by the chiral and (light-quark) continuum extrapolations are con- 
trolled with rooted staggered chiral perturbation theory (rSXPT) [22, 23] applied to heavy- 
light mesons. In Ref. [39], the heavy-light decay constant was calculated to one-loop in 
rSxPT at leading order in the heavy-quark expansion [(I/Mh)"], where Mh is a generic 
heavy-light meson mass. A replica trick is used in rSXPT to take into account the effect of 
the fourth root of the staggered determinant [15, 21]. 

In addition to using the form calculated in Ref. [39], we also use a chiral fit form that 
includes, in the loops, the effects of hyperfine splittings {e.g., M'^ — Mb) and flavor splittings 
{e.g., Mb^—Mb)- These splittings are ~ 100 MeV, and so not much smaller than M^, despite 
the fact that they are formally of order 1/Mh- Since the lightest pseudoscalar meson masses 
in our simulations are ~ 225 MeV, it is not immediately obvious that including the splittings 
is necessary or useful. Their inclusion is motivated, first of all, by the observation of Arndt 
and Lin [79] that finite-volume effects in the one-loop diagrams can be substantially larger 
with the splittings present. This is mainly due to the fact that accidental cancellations in 
finite volume effects between different diagrams at {1/Mh)'^ disappear once splittings are 
included. As described below, it is not difficult to include the splitting effects into the 
calculation of Ref. [39]. We also discuss the extent to which including the splittings, but not 
other effects that could occur at order 1/Mh, is a systematic approximation. In practice, 
we do fits both including and omitting the splittings, and use the difference as one estimate 
of the chiral extrapolation error. For central values, we include the splittings, because this 
yields a more conservative estimate of finite-volume effects. 

With staggered quarks, the (squared) pseudoscalar meson masses are 

Mfj,^^ = Bo{ma + rub) + a^Ag, (6.1) 

where rria and are quark masses, Bq is a parameter of XPT, and the representation of the 
meson under the taste symmetry group is labeled by ^ = P, A, T, V, I [22] . The exact non- 
singlet chiral symmetry of staggered quarks as m^, nib ^ ensures that Ap = 0. All of these 
pseudoscalars appear in the "pion" cloud around the heavy-light meson in the simulation, 
and all of them therefore affect the decay constant. 

Working at leading order [{1/Mh)^] in the heavy-quark expansion and at one loop, or 
next-to-leading order (NLO), in the chiral expansion, the rSXPT expression for the decay 
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constant with light valence quark q takes the form [39] 



E 



d 



(2,q) V,-' 



d 



dM^r 



(6.2) 



where is the light valence-quark mass, e runs over the sea quarks, the lighter two of 
which have masses mi, and the heavier, rrih.^ The parameter (p^ is independent of the light 
masses, and p is an analytic function. We fit the charm and bottom systems separately, so 
05/ depends, in practice, on the heavy-quark mass. The meson mass Mq ,e is similar to Mab,s^ 
in Eq. (6.1), but constructed from a valence quark-antiquark, qq. The light-meson decay 
constant f = 130.4 MeV and the H-H*-tt coupling g^^ controls the size of the one-loop 
effects. Taste-violating hairpin diagrams, which arise only at non-zero lattice spacing, are 
parameterized by 5^ and 5y. The residue functions I^J''^\{M} , {/i}) are defined in Ref. [23]. 
Chiral logarithms are written in terms of the functions £(M^) [80]: 



^{M'^) = M^ln 
i{M^) = 
6i{ML) = 



In — + 5i (ML) j 



[infinite volume], 

[spatial volume L^] 



4 ^Ki{\r\ML) 



(6.3) 
(6.4) 
(6.5) 



Here is the chiral scale, Ki the Bessel function of imaginary argument, and r any non-zero 
three- vector with integer components. The mass sets in the residue functions of Eq. (6.2) 
are 



{Ml,, Ml}, 

■ {m'q,mIm'^,}, 



(6.6) 
(6.7) 

(6.8) 



where Mjj (Ms) is the mass of the pseudoscalar // (hh) meson. 

The salient feature of the chiral extrapolation of (pHq is that the chiral logs have a charac- 
teristic curvature as nig — )■ [81]. At non-zero lattice spacing, the presence of the additive 
splittings a^Aj in the meson masses reduces the curvature of the chiral logarithms. The 
characteristic curvature returns, however, as the continuum limit is approached. 



^ The physical values of the average up-down quark mass and of the strange-quark mass are denoted by 
m = (m„ -f- md)/2 and m^, respectively. 
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To combine data from several lattice spacings into one chiral extrapolation, it is necessary 
to convert lattice units to (some sort of) physical units. As mentioned in Sec. Ill A, we 
convert in two steps, first by canceling lattice units with the appropriate power of ri/a. In 
particular, pseudoscalar meson masses [c/. Eq. (6.1)] become rjM^^^^ = {ri/ aY{aMab,^Y i 

and the decay constant [c/. Eq. (6.2)] becomes r\^'^(j)H = (ri/a)^/^(a^/^0i^), with a^/'^(j)H 
determined from Analyses I or II (c/. Sec. IV). Strictly speaking, one must take the quark 
mass dependence of ri into account, either separately or by modifying the right-hand side of 
Eq. (6.2) accordingly. At the present level of accuracy, we ignore this subtlety, canceling units 
ensemble-by-ensemble with the computed ri/a. Since ri is expected to depend smoothly on 
nil and m/i, we are unlikely to introduce an uncontrolled error into the extrapolated decay 
constants. (After completing the chiral-continuum extrapolation in ri units, we then use 
ri = 0.3117(22) fm (c/. Sec. Ill A) to convert to MeV.) 

To quantify the size of NLO (and higher) corrections to XPT, it is useful to define di- 
mensionless parameters Xq, xi and Xh proportional to the quark masses rrig, mi and nih- 



Xq,l,h 



{riBo){ri/a){2amq^i^h) 



Since the splittings a^Ag are added to the quark mass terms in Eq. (6.1), it is similarly 
useful to define 



^A, - (6.9) 



r^a^A 

Xa = g \ 2 ' (6-10) 



where A is the average pion splitting 

A = ^(Ap + 4A^ + 6AT + 4Ay + A/). (6.11) 

The Xi are in "natural" units for XPT, in the sense that one expects that chiral corrections, 
when written as series in the x,, have coefficients [or low-energy constants (LECs)] that are 
of order 1 . 

We then take the analytic function p in Eq. (6.2) to have the following form at NLO 

^val(Xg + XA^J + /^sea(2X/ + Xh + S^A.^a) + La ^ „ ^ ,^ 4 , (6.12) 

where Lvai, -Z^sea and La are quark-mass-independent LECs that we fit from our data, and 
we define 

9 4 , , 

XA.,, = -^Xa - -XA,, (6.13) 

9 2 

The low-energy constants Lvai, -Z^sea and La depend implicitly on the chiral scale A;^, so that 
the complete expression, Eq. (6.2), is independent of A^. As in Ref. [39], we choose to 
include the dependent terms xa^^^ and xa^ i in the coefficients of L^ai and Lgea so that 
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these coefficients represent those combinations of meson masses that arise naturally under 
a change of in the chiral logarithms. 

The LEG La arises from analytic taste-violating effects; it serves as a counterterm to 
absorb changes proportional to the taste-violating hairpins 5^ and 5y under a change in 
chiral scale. As such, we take the a? coefficient of La in Eq. (6.12) to vary with lattice spacing 
like x/^^^y As long as then appears as an independent ffi parameter, the introduction of 
the XAsea xa^^i terms in the coefficients of L^ai and L^ea. in Eq. (6.12) has a negligible effect 
on the results from the chiral ffis. However, we find that the introduction of these terms 
significantly reduces the magnitude of L^; in other words, most of the discretization error 
from the light quarks appears to be due to the a? dependence of the light meson masses in 
the chiral loops. We leave Lvai , -^^sea and La unconstrained in the ffis that determine central 
values; their size is of 0(1) as expected (and is in fact < 0.6). 

In the region of the strange-quark mass, the data for the decay constants show some 
curvature, and at least some quadratic terms in the quark masses (NNLO effects) must in 
general be added in order to obtain acceptable (p > 0.01) ffis. There are four such LECs, 
giving a NNLO contribution to p of the form 

Qixl + Q2(2a;z + Xhf + Q3Xg{2xi + Xh) + Q4{2xf + xl). (6.15) 

Fits omitting the Qi and Q3 terms give poor confidence levels and are rejected; adding the 
Q2 and Qi terms does not change the ffi results much, but increases over-all errors by up to 
30%. To be conservative, we include all four terms in ffis for central values; other acceptable 
fits (for example, fixing Q2 or Q4 or both to zero) are included among the alternatives used 
to estimate the systematic error of the chiral extrapolation. 

For the central-value ffis, the Qi are mildly constrained by Gaussian priors with central 
value and width 0.5, since that is roughly the expected size in natural units. After ffiting, 
the posterior values satisfy \Qi\ < 0.5, and Qi and have errors ~ 0.05 (much less than 
the prior width), indicating that they are constrained by the data. Q2 and Q^ have errors 
~ 0.5, indicating that they are largely constrained by the priors. Changing the prior widths 
for the Qi to 1.0 has a negligible effect on central values and errors of the decay constants, 
although the posterior Q2 and Q^ typically increase in size and error, as expected. 

While the chiral form introduced so far gives acceptable simultaneous ffis to our data 
from all available lattice spacings, we still need to estimate the size of heavy-quark and 
generic light-quark discretization errors. Following the Bayesian approach advocated in 
Refs. [68, 69], we add constrained lattice-spacing-dependent terms to the ffi function until 
the statistical errors of the results cease to increase appreciably. For the heavy quark, we take 
up to six such terms, fEijnoo), fx{rnoa), fyirnoa), /^(moa), fs{mQa), and /2(^«oa), where 
mo is the heavy quark bare mass. Details about the origin and form of these six functions are 
given in Appendix A. These functions estimate fractional (not absolute) errors, and as such 
are included within the square brackets in Eq. (6.2) (or its equivalent, Eq. (6.20) below). 
The ffist three are O(a^) corrections and are added to the fit with coefficients Zi (aA)^, 
i G {E, X,Y}, where A is a scale characteristic of the heavy-quark expansion, and the Zi 
are parameters with prior value and prior width 1 (for /y) or a/2 (for and fx, since 
they each appear twice in the analysis of Appendix A). The next two terms are O(asa) 
corrections, added with coefficients ZittgaA, i G {B,3}, with Zi taken to have prior value 
and prior width 1 (for /b) or \/2 (for /a, again because it appears twice). The final term 
arises from the propagation to the decay constants of heavy-quark errors in the tuning of 
the heavy-quark hopping parameter, n. It comes in with coefficient Z2 (aA)^, with Z2 having 
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prior value and prior width 1. We take a large value A = 700 MeV, which provides 
conservatively wide priors, especially for the first five terms. Once one of each of the first 
two types of terms is added, the errors already reach ~ 80% of their values with all six 
added. 

Similar terms representing generic light-quark errors, which are not automatically in- 
cluded in the fit function (unlike taste- violating terms), may also be added. With the asqtad 
staggered action, generic discretization effects are of O(asa^). We allow the physical LECs 
05/5 -^vah -^sea, Qi, Q2, Q3, and Q4, to havc small relative variations with lattice spacing 
with coefficients Cias{aA)'^ , where i stands for any of the seven physical LECs, A is again 
taken to be 700 MeV, and the Ci have prior value with prior width 1. This corresponds 
to a maximum of about a 3% difference for a given LEG between the a ~ 0.12 fm and 
the a ~ 0.09 fm ensembles. Once several heavy-quark discretization terms have been intro- 
duced, these light-quark terms further increase the total error of individual decay constants 
by less than 10%. However, the errors on the decay constant ratios fos/ fD+ and fsj fB+ are 
significantly increased by light-quark discretization effects, because the heavy-quark effects 
on the ratios cancel to first approximation. For our central values, we include all six heavy- 
quark and all seven light-quark terms, so the total error from a given fit should estimate 
all (taste-conserving) discretization errors, as well as normal statistical effects. To estimate 
"heavy-quark" and "light-quark" discretization effects separately, we set to zero the light- 
or heavy-quark discretization terms, respectively, and then subtract the statistical errors in 
quadrature. Such separate errors are not relevant to any final results quoted below, but are 
included as separate lines in the error budget for informational purposes. 

As mentioned above, our preferred fit form modifies Eq. (6.2) by including the effects 
of hyperfine and flavor splittings of the heavy-light mesons in one-loop diagrams. We now 
briefly describe how one may adjust the results of Ref. [39] to include these splittings. 
In Eq. (6.2), the contributions proportional to come from diagrams with internal H* 
propagators, namely the left-hand diagrams in Fig. 5 of Ref. [39]. Contributions with no 
factor of come from diagrams with light-meson tadpoles, namely the right-hand diagrams 
in Fig. 5 of Ref. [39] . The latter have no internal heavy-light propagators, so are unaffected 
by any heavy-light splittings. The splittings in the former diagrams depend on whether the 
light-meson line is connected (Fig. 5a, left, of Ref. [39]), or disconnected (Fig. 5b, left). In 
the disconnected case, the H* in the loop always has the same flavor (g) as the external Hq, 
so there is no flavor splitting between the two, only a hyperfine splitting. In the connected 
case, the H* in the loop has the flavor of the virtual sea quark loop (which we labeled by e 
in Eq. (6.2)), so there is flavor splitting with the external Hq, in addition to the hyperfine 
splitting. 

Let A* be the (lowest-order) hyperfine splitting, and S^q be the flavor splitting between 
a heavy-light meson with light quark of flavor e and one of flavor q. At lowest order, is 
proportional to the quark-mass difference (or light-meson squared mass difference), which 
can be written in terms of a parameter Ai: 

5,q ^ 2\,Bo{m, - TTlq) = \^{Ml - M2), (6.16) 

where Me is the mass of an ee light meson. Here we have used the notation of Arndt and 
Lin [79] and included a factor of Bq in the middle expression; Bq is omitted in the notation 
of Ref. [82], Eq. (16), and of Ref. [39], Eq. (45). 

By convention, the mass of the external H is removed in the heavy quark effective theory, 
so the mass shell is at fc = 0, where k is the external three-momentum. When there 
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is no splitting, the internal H* has its pole at the same place, which makes the integrals 
particularly simple, giving the chiral log function i{M^), Eq. (6.3). If a splitting A is present, 
the integrals involve a significantly more complicated function, which we denote 

J(M, A) = (M^ - 2A2) log(MVA^) + 2A'^ - AA'^F{M/A) [infinite volume]. (6.17) 

Here the function F is most simply expressed [83, 84] 



F{l/x) 



tan 



-1 



2 ""^^ 

^Sln(x + y^^^), 



if \x\ < 1, 
if \x\ > 1, 



(6.18) 



which is valid for all x. 

It is then straightforward to write down the generalization of Eq. (6.2) to include split- 
tings. The basic rule is to replace 



i{M^) J(M, A) 



(6.19) 



in the terms proportional to g^. It is not hard to show that J(M, 0) = £(M^), so this 
replacement is consistent with the original result neglecting the splittings. In making the re- 
placements, one must choose the correct value of the splitting A in each term. As mentioned 
above, in terms that come from the diagram with a disconnected light-meson propagator, one 
must put A = A*. But in terms that come from the diagram with a connected light-meson 
propagator, one must put A = A* + Seg, because the internal heavy-light meson is a if*, 
while the external meson is an Hg. The result for the heavy-light meson decay amplitude 
including the splittings is then 



1 + 



1 
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i?f '1(-M?''^; /ii/'V(M„ A*) +[V^A] 



+ p{mg,mi,mh,a^ 



(6.20) 



It is also straightforward to include finite- volume effects into Eq. (6.20). One simply 
replaces 

J{M,A) ^ J{M,A) + 6J{M,A,L), (6.21) 
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where SJ{M, A, L) is the finite-volume correction in a spatial volume L^. The correction 
can be written in terms of functions defined in Refs. [41, 79]: 



2 



M' 

SJ{M,A,L) = 6i{ML) - IGtt^ 

3 



2A A^ - 

—Jfv{M, a, L) + Kpv{M, A, L) 



(6.22) 



with 

KMM,A,L) = ^JMM,A,L), (6.23) 

and 6i{ML) as given in Eq. (6.5). 

Before turning to the fit details and results, we briefiy discuss the extent to which in- 
cluding the splittings as in Eq. (6.20), and not other possible 1/Mh effects, is a systematic 
improvement on Eq. (6.2). In fact, in a parametric sense within the power counting intro- 
duced by Boyd and Grinstein [82], this is a systematic improvement, as long as we make 
some further specifications as to how Eq. (6.20) should be applied. As we detail below, 
however, the power counting of Ref. [82] is only marginally applicable to our data. For that 
reason we ultimately fit to both Eq. (6.20) and Eq. (6.2) and take the difference as one 
measure of the chiral extrapolation error. 

For the following discussion, let A be a generic splitting (A* or 6eq or a linear combination 
of the two), and M be a generic light pseudoscalar mass. The power counting introduced in 
Ref. [82] takes 

A2, AM, M2 ^ ^ ^ 

— '——^ < A ~ M. (6.24) 

For our data, treating A and M as the same size is not dangerous, even though A is 
significantly smaller than our simulation M values — at worst this means that we include 
some terms unnecessarily. The condition M'^/Mh ^ A, which is necessary to drop other 
1/Mh contributions as still higher order, is marginally valid, however. For the D system, 
M'^/Md ^ 130 MeV, which is roughly of the same size as A* and 6sd- For the B system, 
Mj^/Ms ~ 47 MeV, of the same size as A* but somewhat less than Sgd- For the purposes of 
the chiral extrapolation, however, what matters is the applicability of the power counting 
at the lowest simulated light meson masses, not its applicability at M^-"^ For our lightest 
simulated pions with mass ~ Mi^/2, we can reduce the left hand side of the inequality in 
Eq. (6.24) by a factor of four, at which point it becomes reasonably applicable. 

Having tentatively accepted the power counting of Eq. (6.24), it is clear that F{M/A) in 
Eq. (6.17) should be treated as 0(1). Then the difference between J(M, A) and the chiral 
logarithm it replaces, £(M^) is of the same order as ^(M^) itself, so including the splittings 
becomes mandatory at the one-loop order to which we are working. The next question is 
whether Eq. (6.2) includes all effects to this order. As discussed by Boyd and Grinstein, the 
key issue is whether operators with two or more derivatives (two or more powers of residual 
momentum k) on the heavy fields can contribute. Such operators are suppressed by 1/Mh 
relative to the leading-order heavy-light Lagrangian, which has a single derivative. Since 
we are keeping A*, which is also in principle a 1/Mh effect, one might worry that such 
operators could contribute at the same order. The power counting implies, however, that 
the relevant diagrams pick up a factor of {A, M)/Mh relative to the terms being kept in 



^ We assume here that the fit to the data is good over the full mass range simulated. It is not important 
for the chiral extrapolation that the fit be systematic in the region around il/^, but it must describe the 
data in that range so that we can correctly interpolate to the physical kaon mass. 
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Eq. (6.20). The reason for the difference is that the exphcit extra factor of k turns into A or 
M — the only dimensional constants available — after integration. In the term that generates 
the hyperfine splitting itself, in contrast, the dimensional quantity balanced against 1/Mh is 
A — a heavy-quark QCD scale — rather than M. The power counting in Eq. (6.24) effectively 
treats A as larger than M (so that A ~ h? /Mu ~ M). Similarly, the term that generates 
the flavor splittings has a single factor of and no residual momentum, and Eq. (6.24) 
effectively takes niq ~ fc in such terms. 

Boyd and Grinstein do find some other contributions at the same order as Eq. (6.20), but 
most come from terms that are simply KjMu times terms in the leading-order heavy-light 
Lagrangian or current, and thus give simply an overall factor times the result without them. 
The exceptions are the terms multiplied by in Eq. (15) of Ref. [82] and by pi in Eq. (18) 
of Ref. [82] . These are operators that have the same dimension as the original Lagrangian 
current operators, but that violate heavy-quark spin symmetry, and therefore give different 
contributions to the pseudoscalar and vector meson decay constants at this order. Since 
we are only looking at pseudoscalar meson decay constants here, however, and since these 
effects are flavor-independent, we can also absorb all of the 1/Mh effects into (1) the effects 
of the splittings in the loop, described by Eq. (6.20), and (2) an overall factor in front of the 
full one-loop result. 

The overall factor in Eq. (6.20) is l/(167r^/^). Since / is not fixed at one loop, one should 
in any case allow it to vary over a reasonable range, which we take to be to fx- We allow 
such variations even when we fit to the form without splittings, Eq. (6.2). The difference 
between using and fx corresponds to a 45% change in the size of the one-loop coefficient, 
but produces only a 1 to 3 MeV change in the decay constants."^ We therefore assume that 
any further 1/Mh uncertainty in l/(167r^/^) has negligible effects on our results. 

Finally, there is a question of whether terms coming from taste violations contribute 
something new at the same order in which we include splittings. Since taste-violating terms 
in the Lagrangian can enter just like light-quark masses, this is a possibility in principle. 
Corresponding to the terms in the quark masses that generate flavor splittings of heavy- 
light mesons (c/. Eq. (45) of Ref. [39]), there are taste-violating terms given in Eq. (51) 
of that paper. Just as for the quark-mass terms, however, we are only interested here in 
contributions that change the heavy-light meson mass, not ones coupling the mesons to pion 
fields. When the pion fields are set to zero, all the terms in Eq. (51) of Ref. [39] just give a 
constant heavy-light meson mass term proportional to a? that contributes equally to the H 
and H* masses of all valence flavors. Terms that produce a hyperfine splitting would have 
to also violate heavy quark spin symmetry, and hence be of order o^K/Mh- Similarly, terms 
that produce flavor splitting would need to violate flavor symmetry, and hence be of order 
a^mq/A^. Both such contributions are higher order in our power-counting. Since there is 
no splitting, there is no contribution to the decay constants because the effect will vanish 
when we put the external B or D meson on mass shell. 

In our chiral fits, we take the physical light-quark masses, as well as the parameters Bq, 
Ag, and 6y, from the MILC Collaboration's results of rSXPT fits to light pseudoscalars 
masses and decay constants [12, 85] on ensembles that include lattice spacing a ^ 0.15 fm 
through a ~ 0.06 fm. Table XII shows the values used. In general, we use older MILC 
determinations since newer versions, e.g., those in Ref. [66], do not cover the full range of 
lattice spacings employed here (but are consistent where they overlap). The exceptions are 

^ Most of the change in the size of the overaU coefficient is compensated by a change in the LECs that come 
from the fit to our data. 
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TABLE XII. Inputs to our heavy-light chiral fits taken from the MILC Collaboration's light-meson 
chiral fits [12, 85]. The physical bare-quark masses m„, m^, m = (m„ + mrf)/2, and are 
determined by demanding that the charged pion and kaons take their physical masses after the 
removal of electromagnetic effects. Errors in the masses are due to statistics, chiral extrapolation 
systematics, scale determination, and (for and m^) the estimate of electromagnetic effects, 
respectively. "Continuum" values are found from chiral fits that have been extrapolated to the 
continuum, but masses are still in units of the "fine" (a ~ 0.09 fm) lattice spacing, and with the 
fine-lattice value of the mass renormalization. Values for rfa^d'j^ and r^a^S'y take into account newer 
MILC analyses [66] as noted in the text. The light-meson analysis determining these quantities 
assumes that they scale like the taste- violating splittings and are larger by a factor of 1.68 on 
the 0.15 fm lattices than on the 0.12 fm lattices, and smaller by a factor 0.35 on the 0.09 fm lattices 
than on the 0.12 fm lattices. The statistical and systematic errors on tiBq and r\o?lS.^ are not 
given here; such errors have negligible effect on the heavy-light decay constants. 



Quantity 




Lattice 


spacing 






a w 0.15 fm 


a ^ 0.12 fm 


a ^ 0.09 fm 


"continuum" 


arus X 10^ 


4.29(1)(8)(6) 


3.46(1)(10)(5) 


2.53(0)(6)(4) 


2.72(1)(7)(4) 


am X 10^ 


1.55(0)(3)(2) 


1.25(0)(4)(2) 


0.927(2)(27)(13) 


0.997(2)(32)(14) 


amd X 10'^ 


2.20(0)(4)(3)(5) 


1.78(0)(6)(3)(4) 


1.31(0)(4)(2)(3) 


1.40(0)(5)(2)(3) 


aruu X 10^ 


8.96(2)(17)(13)(49) 


7.31(2)(23)(10)(40) 


5.47(1)(16)(8)(30) 


5.90(1)(19)(9)(32) 




6.43 


6.23 


6.38 


6.29 


rfa^A^ 


0.351 


0.205 


0.0706 







0.555 


0.327 


0.115 







0.721 


0.439 


0.152 





rfa^Ai 


0.897 


0.537 


0.206 









-0.28(6) 











0.00(7) 








the values of the taste- violating hairpin parameters rfa'^S'^^ and rfa'^Sy. For them, the newer 
analysis including two-loop chiral logarithms gives larger systematic errors and a changed 
sign of the central value of rfa'^6y, which has always been consistent with zero. For these 
parameters, we therefore use the wider ranges listed in Table XII, which encompasses both 
types of analyses. For comparison, the results of the analysis of Ref. [85] were rfa^5^ = 
-0.30(1)(4) and rfa^d'y = -0.05(2)(4). 

In order to fit Eq. (6.20) to our lattice data, it is also necessary to input values for the 
hyperfine splitting A* and for Ai in Eq. (6.16). For B mesons, we have [51] 



A* = Mb* -Mb^ 45.8 MeV, (6.25) 
= Mb, -Mb^ 87.0 MeV, (6.26) 
Ai ^ 0.192 GeV~\ (6.27) 



where we use Me = Ms = 0.6858(40) GeV [67] and Mq = M^o ^ 135.0 MeV to obtain Ai 
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from the experimental data. Similarly, for D mesons, we have 

A* = Md* - Mdo ~ 142.1 MeV, (6.28) 
Ssd = Md^ - Md^ ~ 98.9 MeV, (6.29) 
Ai ^ 0.219 GeV'K (6.30) 

In the chiral fit, we input the relevant physical A* and Ai from either Eqs. (6.25)-(6.27) 
or (6.28)-(6.30), and then use Eq. (6.16) with the actual mg and rrig from each data point, 
and Bo the slope for a given ensemble, from Table XII. We emphasize here that Bq comes 
from a simple tree-level chiral fit of light meson masses to Eq. (6.1). This is adequate for 
our purposes, since the resulting meson masses are only used within the one-loop chiral 
logarithms. 

We can now present the actual chiral fits and show how we extract results and systematic 
errors from them. Recall that we compute (pHq for many combinations of the valence and light 
sea-quark masses, and at three lattice spacings: a ~ 0.15, 0.12, and 0.09 fm. We fit all the 
decay constant data to the form given either by Eq. (6.20) or by Eq. (6.2). One-loop finite- 
volume effects are included through Eq. (6.21) or Eq. (6.4). There are four unconstrained 
free parameters in our fits: the LO parameter 0^, and the one-loop LECs Lvai, -Z^sea, La 
[Eq. (6.12)]. The central fit fixes the chiral coupling / at f^, but a range of couplings are 
considered in alternative fits, as described in more detail in Sec. VII. Similarly, the H-H*-7r 
coupling Qtt, which is poorly constrained by our data, is taken in the range 0.51 ± 0.20. 
This encompasses a range of phenomenological and lattice determinations [83, 86-91], as 
discussed in Ref. [53]. In the central fit, g.,^ is held fixed at 0.51, while it is varied in 
alternative fits described in Sec. VII. Although changing g^^ is equivalent to changing / 
when splittings are omitted [cf. Eq. (6.2)], the effects are inequivalent when splittings are 
included [cf. Eq. (6.20)]. This is especially true of the finite-volume effects, for which the 
splittings have the potential to produce significant changes [79]. 

Some additional parameters constrained by Bayesian priors are also included in the chiral 
fits, as discussed above. The taste-violating hairpin parameters 6y and 5^ are given by 
the ranges in Table XII. In addition, up to six heavy-quark and up to seven light-quark 
lattice-spacing dependent terms, are added for investigation of discretization effects. Except 
where otherwise noted, all twelve such terms are included in the fits plotted below: this 
gives errors that include true statistical errors plus our estimate of discretization effects 
from the heavy quarks and generic (taste non-violating) discretization errors from the light 
quarks. In addition, some or all of the (mildly) constrained NNLO LECs, Qi, . . . ,Q4,, are 
included. Again, unless otherwise noted, the fits below include all four such parameters; 
such fits tend to give larger (and hence more conservative) errors than fits that restrict the 
number of these parameters. In total, there are 23 fit parameters in the central fits: the 19 
constrained parameters listed in this paragraph, and the 4 unconstrained parameters listed 
in the previous paragraph. 

B. Chiral fits and extrapolations for the D system 

Figure 3 shows our central chiral fit to r^^^0D+ and r^^"^ (pr,^. Data from ensembles at 
a 0.15 fm, a ~ 0.12 fm and a ~ 0.09 fm are shown, but the a ~ 0.15 fm ensembles 
are not included in the fit. The points and co variance matrix are obtained from Analysis I 
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(Sec. IV B) of the two-point functions. For clarity, only the unitary (full QCD) points are 
shown for 0^) (and approximately unitary for ^Dj,), but the fit is to all the partially-quenched 
data on the a ~ 0.12 fm and a ~ 0.09 fm ensembles. The fit properly takes into account the 
covariance of the data; x^/dof and the p value (goodness of fit) are reasonable, as shown. 
The points in Fig. 3 are plotted as a function of mass m^,, where, for ipD+i the light valence 
mass TUq and the light sea mass mi are given by = mi = mx. For only mi = mx 
varies, while is held fixed at the value m^^ near the physical strange mass m^.^ In order 
to be able to compare ensembles at different lattice spacings, we have adjusted the bare 
quark masses by the ratio Zm/ Z'^^^'^, where Zm is the (one-loop) mass renormalization 
constant [92], and Z^*^^^™ is its value on the a ~ 0.09 fm ensembles. 

The continuum extrapolation is carried out by taking the fitted parameters and setting 
= in all taste- violating terms (parameterized by A^, 5^, 5y, and L^), all heavy-quark 
discretization effects (parameterized by ze^ zx, zy, zb, z^, and Z2) and all generic light-quark 
discretization effects (parameterized by Ci). The red lines (solid for 4>d+, dotted for 0/5 J 




0.00 0.05 0.10 

^ ^ f'V /r/O.OQfmx 

^1 ^ (2m/Zm ) 

FIG. 3. Central chiral fit for the D system, based on Analysis I of the fits to 2-point correlators. 
Only (approximately) unitary points are shown. Data from ensembles at a w 0.15 fm, a ~ 0.12 fm 
and a ~ 0.09 fm are shown, but the a ~ 0.15 fm ensembles are not included in the fit. The bursts 
show extrapolated values for 0/)^ and (po+y with the purely statistical errors in bright red and the 
statistical plus discretization errors in darker red. The physical strange-quark mass corresponds to 
an abscissa value of nix ~ 0.1. 

^ On the a « 0.15 fm ensembles, rns„ is equal to the value of the strange sea quark mass ruh {airis^ = 
0.0484), but on the other two ensembles we take it lower than m/i, because rrih has been chosen somewhat 
larger than the physical strange mass. In the figure, arUs^ = 0.415 for the a ~ 0.12 fm ensembles and 
ams„ = 0.272 for the a « 0.09 fm ensembles. 
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show the effect of extrapolating to the continuum and setting the strange quark mass (both 
sea, m/i, and valence, m^^) to the physical value m^. 

Finally, the bursts give the result after the chiral extrapolation in the continuum, i.e., 
setting = rrid for (j)D+, and rrix = rh for (p^^. The larger, dark red, error bars on the 
bursts show the total error from the fit, which includes heavy-quark and generic light-quark 
discretization errors using Bayesian priors, as described above. The smaller, bright red error 
bars, show purely statistical errors, which are computed by a fit with all the discretization 
prior functions turned off. In plotting the red line for (j)D+, the light sea mass is shifted 
slightly {nil = f^x + rh — rrid) so that it takes its proper mass when rrix = rrid. (We neglect 
isospin violations in the sea.) The small mass differences between m and rrid (and the 
corresponding difference between rh and rriu for the B~^) produce changes in that are 
much smaller than our current errors, but we include them here with an eye to future work, 
where the precision will improve. 

The trend of the data for the coarsest lattice spacing (a ^ 0.15 fm, the magenta points 
in Fig. 3) tends to be rather different than for the finer lattice spacings, especially for the 
Dg, which is why we exclude the a ^ 0.15 fm data from the central fit. This trend is even 
more exaggerated for the B system, but with particularly large statistical errors; see Fig. 6 
below. Nevertheless, the effect of including the a ^ 0.15 fm points in the fit is a rough 
indication of the size of discretization errors. Figure 4 shows what happens to the fit when 
these points are included: 0d+ and (pDs each move up an amount comparable to (but less 
than) the size of the larger (dark red) error bars, which represent heavy and generic light 
quark discretization errors (as well as statistical errors, which are smaller). The consistency 
is reassuring. 

As discussed in Sec. IV, we also examine Analysis II of the 2-point functions. Figure 5 
shows the effect of using Analysis II in the chiral fits. The differences in the decay constant 
results between Fig. 3 and Fig. 5 are included in the decay-constant error budgets as a 
"fitting error". Note that the covariance matrix calculation in Analysis II results in an 
apparent underestimate of (and, consequently, a high apparent p value). We believe 
that this stems from binning of the data to remove autocorrelation effects, which has the 
disadvantage of reducing the number of samples used to compute the covariance matrix. 
It is then difficult to determine small eigenvalues accurately. Indeed the eigenvalues of the 
(normalized) correlation matrix tend to have a lower bound of ~ 10~^ to 10"^ with this 
approach, whereas they typically go down to 10~^ in Analysis I. [Recall that in Analysis I 
we keep all samples, and deal with autocorrelation effects by Eq. (4.13).] Nevertheless, the 
difficulty with small eigenvalues explains only a small fraction of the difference between the 
results from Analyses I and II. For example, fo is changed by only 0.2 MeV when we smooth 
eigenvalues from Analysis I that are less than 10~^, following the method of Ref. [93]. This 
may be compared to the total difference between fo in Analyses I and II, which is 1.7 MeV. 

C. Chiral fits and extrapolations for the B system 

Results for the B system closely resemble those for the D system in most respects. One 
important difference is that the signal-to-noise ratio is worse for the B system because the 
mass difference that controls the noise, 2m b — rrirj^ — m^, increases with the mass of the 
heavy quark [94]. Therefore, the preferred fit in Analysis I for the charm case (1 simple 
exponential -|- 1 oscillating exponential at large tmin) is too noisy here, and we must use fits 
with an extra excited state and smaller tmin (see Sec. IV B). Consequently, our S-system 
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results have larger statistical errors. On the other hand, heavy-quark discretization errors 
are smaller in the B system. In the HQET analysis of discretization effects they appear in 
the heavy-quark expansion, which works better for 5's to begin with [60]. 

Figure 6 shows, for unitary points only, our central chiral fit for the B system. This 
is based on Analysis I of the 2-point functions. As in Fig. 3, the red lines (solid for 
dotted for 0bJ show the effect of extrapolation to the continuum and setting the strange 
quark mass to its physical value m^. For the solid red line, the light sea mass is again shifted 
slightly, but now mi = rux + rh — ruu, so that it takes its proper mass when rrix = rriu. The 
bursts show the final results, and come from setting rUx = triu for 0^+ and rrix = rh for 0^^. 
As before, the smaller, bright red, error bars, show purely statistical errors, and the larger, 
dark red, error bars come from the fit with Bayesian priors and include heavy-quark and 
generic light-quark discretization errors as well as statistical errors. 

In Fig. 6, the a ~ 0.15 fm data are both noisy and far from those of the finer lattice 
spacings. Therefore, these ensembles are again dropped from the central fit. Figure 7 
shows the effect of including the a ~ 0.15 fm points. Note that the resulting continuum- 
extrapolated line for (dotted red line) now has what appears to be a rather unphysical 
shape, showing a significant initial increase as the light sea-quark mass is decreased, starting 
at the right side of the graph. Hence, the differences caused by including the a ^ 0.15 fm 
points is 10 to 20% larger than the dark red error bars in Fig. 6, and 40 to 60% larger than 
discretization errors estimated by removing the statistical errors from the dark red bars. 
Because the trend for a ~ 0.15 fm is so different from the other spacings, and because of 
the unphysical behavior when these points are included in the fit, we believe this difference 




« exlrap J gtal. err. J w/ di^cretiz. err. 



0.00 0.05 0.10 

X (Z^/Zr^-) 

FIG. 4. Same as Fig. 3, but including points at a ~ 0.15 fm in the chiral-continuum fit. 
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overestimates the true discretization error, and we do not enlarge the errors coming from 
the fit. 

Figure 8 shows the effect of using Analysis II of the correlation functions. In order to 
make these comparisons as direct as possible, we first turn off all the Bayesian discretization 
terms in the fits. Compared to the results from Fig. 6, this fit gives a value of fs^ about 1 
MeV higher and a value of /b+ about 2 MeV lower. These differences are included in our 
estimate of the fitting errors due to excited state contamination in Sec. VII. 

VII. ESTIMATION OF SYSTEMATIC ERRORS 

In this section, we present a careful, quantitative accounting for the uncertainties in our 
calculation. We consider in turn discretization errors, fitting errors, errors from inputs ri 
and quark-mass tuning, renormalization, and finite-volume effects. Table XIII details our 
error budget. 

A. Heavy-quark and generic light-quark discretization effects 

As described in Sec. VI and Appendix A, we parameterize possible heavy-quark and 
generic light-quark discretization effects and follow a Bayesian approach in including such 
effects in our chiral fitting function. Consequently, the raw "statistical" error that comes 
from our fits is not a pure statistical error but includes an estimate of the errors coming 

0.8 
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X (Z^/Zr^-) 

FIG. 5. Same as Fig. 3, but using Analysis II of the 2-point function. 
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D; 6.12fm/0.09fm fi't 




« exlrap J gtal. err. J w/ di^cretiz. err. 



from the discretization effects. This inclusive error is shown with the dark red error bars in 
the plots in Sec. VI, and is listed in the first line of Table XIII. 

For informational purposes, it is useful to break down this inclusive error into its compo- 
nent parts, at least approximately. We can see what errors to expect and, hence, target for 
improvement in future simulations. In particular, with our current actions, the light-quark 
and heavy-quark discretization errors should behave differently as a function of lattice spac- 
ing, with heavy-quark errors decreasing more slowly as a is reduced. To extract the pure 
statistical errors, we rerun the fits with all the Bayesian discretization terms set to zero. 
We then find the pure heavy-quark (or pure light-quark) discretization contributions, by 
turning back on the heavy-quark (light-quark) terms, and then subtracting in quadrature 
the pure statistical errors from the resulting raw errors. These individual errors are shown 
in Table XIII in parentheses. Note that the total error at the bottom of the table includes 
the error on the first line, not the sum of the three errors in parentheses, when these differ. 
Note also that the discretization errors are similar to what we would have obtained with less 
sophisticated power counting. 
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FIG. 6. Central chiral fit for the B system, with data from Analysis I of the 2-point functions. 
Only (approximately) unitary points are shown. Data from ensembles at a ~ 0.15 fm, a ~ 0.12 fm, 
and a ~ 0.09 fm are shown, but the a ~ 0.15 fm ensembles are not included in the fit. The bursts 
show extrapolated values for (f)B^ and </)s+, with the purely statistical errors in bright red and the 
statistical plus discretization errors in darker red. The physical strange-quark mass corresponds to 
an abscissa value of m-r w 0.1. 
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B. Chiral extrapolation and taste- violating light-quark discretization effects 

As described in Sec. VI, we modify the chiral fit function in a variety of ways to estimate 
the error associated with the chiral extrapolation: 

Xl. Set the chiral coupling / to fx instead of 

x2. Allow the chiral coupling / to be a Bayesian fit parameter, with prior value and 
prior width equal to fK — f-n- 

x3. Replace the H-H*-tt coupling (7^ (which is 0.51 in the central fit) with 0.31 or 0.71, 
which are the extremes of the range discussed in Sec. VI. 

x4. Allow to be a constrained fit parameter, with prior value 0.51 and prior width 0.20. 

x5. Fix to zero those NNLO analytic terms [Q2 and/or Q4 in Eq. (6.15)] that may be 
eliminated without making the fit unacceptably poor. 

x6. Use the chiral function without hyperfine and flavor splittings, i.e., use Eq. (6.2) 
instead of Eq. (6.20). 

x7. Use combinations of modifications xl and x3 or modifications x2 and x3. These 
choices can produce significantly larger deviations since changes in g^^ have a similar 
effect on the fit function as changes in /. 
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FIG. 7. Same as Fig. 6, but including points at a ~ 0.15 fm in the fit. 
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These modifications typically change the decay constant by 1-3 MeV, and the ratios by 
1-1.5%. We take the chiral extrapolation error of a given quantity to be the largest change 
(of either sign) under the above modifications, and list it in Table XIII. In several cases, 
(/d+) IdJ Id+i and f'Bsl fB+) the largest change comes from modification x6, eliminating 
the heavy-light splittings. The fit without the splittings is shown for the D system in Fig. 9. 
It may be compared to Fig. 3 to see the effects: the curvature at small mass for is 
slightly greater without the splittings, which results in a decrease of /^n- of 3.2 MeV. Note 
that the p values of the two fits are almost identical, so the goodness-of-fit cannot be used 
to choose one version of the chiral extrapolation over the other. 

Modifications of / and/or g^^ produce the largest changes in the other quantities, namely 
fosi fB+ and Jbs- In particular, putting f = fx and g.,, = 0.31 results in an increase of +2.9 
for fB+ and +2.8 MeV for Jb^- The modified fit is shown in Fig. 10, and may be compared 
with Fig. 6 to see the effects of the changes. Increasing / and decreasing g.,, both suppress 
the chiral logarithms [cf. Eq. (6.20)] and give fit functions with less curvature and smaller 
slope at low quark mass. 

Since the rSXPT fit functions in Eqs. (6.2) and (6.20) explicitly include one-loop dis- 
cretization effects coming from taste violations in the (rooted) staggered light quark action, 
the chiral error estimates we describe here inherently include taste-violating discretization 
errors. However, it seems unlikely that the current data can accurately distinguish between 
such taste-violating errors of order a^a^ and generic light-quark discretization effects of or- 
der agO^, or even heavy-quark discretization effects. Indeed, the taste-violating LEG La [cf. 
Eq. (6.12)] is not well constrained by our fits and is consistent with zero within large errors. 
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FIG. 8. Same as Fig. 6, but using Analysis II of the 2-point functions. 
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TABLE XIII. Total error budget for the heavy-light decay constants. Uncertainties are in MeV 
for decay constants. The total combines errors in quadrature. The first row includes statistics, 
heavy-quark discretization errors, and generic light-quark discretization errors, as explained in the 
text. Errors in parentheses are approximate sub-parts of errors that are computed in combination. 



Source fD+ (MeV) /d, (MeV) /b+ (MeV) /b, (MeV) 



Statistics © discretization 


9.2 


8.9 


0.014 


5.5 


5.1 


0.013 


(statistics) 


(2.3) 


(2.3) 


(0.005) 


(3.6) 


(3.4) 


(0.010) 


(heavy-quark disc.) 


(8.2) 


(8.3) 


(0.007) 


(3.7) 


(3.8) 


(0.004) 


(light-quark disc.) 


(2.9) 


(1.5) 


(0.012) 


(2.5) 


(2.1) 


(0.007) 


Chiral extrapolation 


3.2 


2.2 


0.014 


2.9 


2.8 


0.014 


Two-point functions 


3.3 


1.6 


0.013 


3.0 


4.1 


0.015 


Scale (ri) 


1.0 


1.0 


0.001 


1.0 


1.4 


0.001 


Light quark masses 


0.3 


1.4 


0.005 


0.1 


1.3 


0.006 


Heavy quark tuning 


2.8 


2.8 


0.003 


3.9 


3.9 


0.005 


uq adjustment 


1.8 


2.0 


0.001 


2.5 


2.8 


0.001 


Finite volume 


0.6 


0.0 


0.003 


0.5 


0.1 


0.003 


Zi/4 and Zt/4 


2.8 


3.4 


0.000 


2.6 


3.1 


0.000 


Higher-order 


1.5 


1.8 


0.001 


1.4 


1.7 


0.001 


Total Error 


11.3 


10.8 


0.025 


8.9 


9.5 


0.026 


The central fits give 
















La 


= +0.6 ±6.5 


{D syst< 


em), 




(7.1) 




La 


= -1.9 ±8.8 


{B systi 


gm), 




(7.2) 



where the error is the raw statistical error. (Note that we do not constrain La by any prior 
width.) The errors in La decrease by about 10% if Bayesian parameters for generic light- 
quark errors are removed, and an additional 10% if the parameters for heavy-quark errors 
are removed. Thus, there is "cross talk" between various error sources, making it difficult 
to completely distinguish the various types of discretization errors. Future work, with more 
and finer lattice spacings, should make a cleaner separation possible. 

C. Fitting errors 

The "fitting errors" are the errors introduced in the analysis of the two-point correlators. 
They represent the effects of various choices of fit ranges and fitting functions, and are 
an estimate of the systematic effect of the contamination by excited states. We compare 
results from the three choices of two-point fitting (see Sec. IV): Analysis I, Analysis II, and a 
modified Analysis I using 1 simple + 1 oscillating state, but values of tmm larger than those 
described in Sec. IV B. 

Some of these differences may, in fact, be due simply to statistical effects, and hence 
already included in the statistical error. Figure 11 shows the differences between values of 
(pBg in Analyses I and II, divided by the average statistical error for each of the common 
partially quenched data points. Only 10 of 74 differences are greater than 1 statistical a. 
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Nevertheless, there appears to be some significant systematic trend in that 46 of 74 points 
are positive. To be conservative, we take the largest difference between the Analysis-I fits 
and the other two fits as the fitting error for each physical quantity, and list it in Table XIII. 
For fi)^ and fs^, the difference is largest for chiral fits based on 2-point Analysis II, while, 
for the other four quantities, the difference is largest for the modified Analysis I. 



D. Scale uncertainty 

We use the scale ri = 0.3117(22) fm to tune the values of the quark masses and convert 
the decay constants into physical units (see Sec. Ill A). To find the scale errors on the final 
results, we shift ri to 0.3139 fm or 0.3095 fm and redo the analysis. Although (pn scales 
like r^^'^^, the change in the results under a change in ri is smaller than pure dimensional 
analysis would imply, because our estimates of the physical light masses and the heavy-quark 

and Kb also shift, producing partially compensating changes in (pH- At ri = 0.3139 fm, 
we shift the light masses in Table XII upward by the scale error shown in that table. [The 
lattice light-quark masses scale like rf, because they are approximately linear in the squared 
meson masses {rirriT^Y and (rim^)^.] Similarly, we shift the tuned Kc and Kb downward by 
the scale error in Table V because the bare heavy quark mass increases with ri. We then 
adjust and at each lattice spacing using Eq. (3.5) and the values of dcp/dn given 
in Table V. Redoing the preferred chiral fits shown in Figs. 3 and 6, extrapolating to the 
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FIG. 9. Same as Fig. 3 but omitting heavy- light hyperfine and flavor splittings in the chiral fit 
function. 
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continuum, and plugging in the adjusted continuum light quark masses gives the scale error 
listed in Table XIII. 



E. Light-quark mass determinations 

To estimate the error from the light-quark mass determination, we follow a similar proce- 
dure to that in the scale-error case. We shift the continuum light quark masses in Table XII 
by the sum in quadrature of all errors except scale errors. This includes the statistical er- 
rors, the chiral errors and, where relevant, the electromagnetic errors. We then plug the new 
masses into the continuum-extrapolated chiral fits and take the difference from the central 
results to give the errors listed in Table XIII. The relative direction of shifts on different 
masses makes little difference in the size of the errors on the decay constants fa^ , fD+ , fs^ , 
and since they are sensitive primarily to the valence quark masses. However, it does 
affect the error of the ratios /ds//d+ and fss/ fB+- The largest effect clearly occurs when the 
strange mass is shifted in the opposite direction from the lighter masses. To be conservative, 
we take the size of change of the ratios in this case as the error, but this is almost certainly 
an overestimate because the statistical and chiral extrapolation errors on the light quark 
masses are positively correlated between the strange mass and the other masses. 

Note that the errors from the light-quark masses in Table XIII are much larger for fo^ 
and than for fD+ and That simply reflects the facts that the decay constants have a 
nonzero limit when the quark masses vanish, and that the dependence on the quark masses 
is reasonably linear. Thus a given percent error in the strange mass produces a much larger 
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percent difference in /^^ and fs,, than the same percent error in the d or u mass does in 
fD+ and fB+. 



F. Bottom and charm quark mass determinations 

The propagation of statistical errors in the tuned and Kb to the decay constants is 
comphcated by the fact that the independent errors at each lattice spacing affect the final 
results in a nontrivial way through the continuum and chiral extrapolations. At each lattice 
spacing, we choose 200 gaussian-distributed ensembles of trial n values with central value 
equal to the tuned values, and standard deviation equal to the statistical error, taken from 
Table V. For a given choice of trial k values at each lattice spacing, we produce an adjusted 
trial data sample by shifting the (pn values according to Eq. (3.5), but with the trial values 
replacing the tuned values. We then perform the complete chiral fit and extrapolation 
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FIG. 11. Difference of (pBq values from Analyses I and II, divided by the average statistical error 
at each of the common valence and sea mass points. The order along the abscissa is arbitrary. 
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procedure on each of the 200 trial data sets. The standard deviation over trials of a given 
decay constant or decay constant ratio is taken to be the heavy quark tuning error, and is 
listed in Table XIII. 



G. Tadpole factor (uq) adjustment 

In order to improve the convergence of lattice perturbation theory, we use tadpole- 
improved actions for the gluons, hght quarks, and heavy quarks [76] . For the gluon and sea- 
quark actions we take the tadpole factor uq from the average plaquette. On the a ~ 0.15 fm 
and a ~ 0.09 fm lattices we use the same choice for the light valence and heavy-quark ac- 
tions. On the a ^ 0.12 fm lattices, however, we use the tadpole factor uq taken from the 
Landau link in the valence-quark action and in the clover term in the heavy-quark action. 
This results in a slight mismatch between the light valence and sea-quark actions on these 
ensembles, and also affects the values obtained for the tuned bottom- and charm-quark 
masses kj, and Kc- The difference between uq obtained from the average plaquette and the 
Landau link is approximately 3-4% on the a ~ 0.12 fm ensembles. 

We propagate this difference through the chiral /continuum extrapolation as follows. 
First, we compute the heavy-strange meson decay amplitudes (pBs and (pOs with both choices 
for Mo on the ensemble with ami/arrih = 0.01/0.05, a ^ 0.12 fm. For each choice of uq, we 
compute (pBs and (pOs directly at the tuned values of Kb and Kc, thereby avoiding an inter- 
polation in K. Next, we renormalize the lattice decay amplitudes using the nonperturbative, 
flavor-diagonal current renormalization factors Zy4 and Zv4 obtained for each case. (We 
neglect the slight difference in the perturbative correction p^^ .) Then, we calculate the 
ratio of the renormalized decay amplitudes, finding no difference within errors: 

^plaquette/^Landau _ 1,005(13), (7.3) 
^plaquette^^Landau _ 1,014(20). (7.4) 

As expected, the Uq dependence from the bare current and renormalization factors mostly 
cancels. Finally, we repeat the chiral/continuum extrapolation shifting 0c and (pb on the 
a ~ 0.12 fm ensembles by the statistical errors reported in Eqs. (7.3)-(7.4). We find that 
these percent-level errors in 0c and (pb lead to approximately 1% errors in the extrapolated 
decay constants and approximately 0.1% errors in the decay-constant ratios. These errors 
are listed as "mq adjustment" in the error budget in Table XIII. 



H. Heavy-light current renormalization 

There are two sources of systematic error in our heavy-light current renormalization. The 
first is due to the perturbative calculation of p^^ and the second is due to the nonpertur- 
bative calculation of Zv-i and Zyi ■ 

The perturbative calculation of p^^ has been carried out to one-loop order. Since p^4^ is 
defined from a ratio of renormalization factors [see Eq. (2.5)], its perturbative corrections are 
small by construction. Indeed, as can be seen from the results for p^4 given in Table VIII, 

we observe very small corrections. For bottom they range from 0.3% at a ^ 0.09 fm to 0.8% 
at a ~ 0.12 fm and 2.8% at a ^ 0.15 fm. For charm they range from less than 0.08% at 
a ^ 0.09 fm to 0.4% at a ^ 0.12 fm and 0.6% at a ^ 0.15 fm. As shown in Ref. [43] the 
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perturbative corrections to the p-factors for the spatial currents, while still small, tend to 
be bigger than those for the temporal currents and V^. We therefore estimate the error 
due to neglecting higher order terms as p|h a^. We take at a ~ 0.09 fm and ~ 0.1, 
which is the largest one-loop coefficient for pyi in the mass range mga < 3. This procedure 
yields a systematic error of 0.7%, which we take for both charm and bottom decay constants. 

The decay constant ratios fsj fB+ and fosl fD+ depend on the corresponding ratios of 
p^4 /p^4 . These ratios differ from unity only because of the small variation of the p^4 
with light valence mass, which is described in Sec. V. We take the variation of the p^ with 

Qq 

light valence mass at a ^ 0.09 fm as the error. This yields an error of 0.1% for both bottom 
and charm. 

The dominant corrections in the heavy-light renormalization factor as defined in Eq. (2.5) 
are due to Zy4^ and which are calculated nonperturbatively. The values (and errors) 
for Zy4 and Zyi are listed in Tables IX and XI, respectively. To obtain the error in 



Zt/4 = . / Zyi Zv4 we add the statistical errors in Zyt and Zyi in quadrature. The 

Qq Y 11 QQ 11 QQ 

error on Zi/4 is dominated by the error on Zyi . The errors are largest, 1.3%, on the 

^Qq ^11 e> ■} ■} 

a ~ 0.09 fm ensemble and they are about the same for both charm and bottom on the two 
finest ensembles used to obtain our main decay constant results. Hence we use 1.3% as our 
estimate for the uncertainty in Zy4 . 

Qq 



I. Finite volume effects 



To study finite volume effects, we use the chiral fit function with heavy-light hyperfine 
and flavor splittings included (Eq. (6.20)), since the effects are known to be larger with the 
splittings than without [79]. The central flt includes the (one-loop) flnite volume corrections, 
Eq. (6.21), on the lattice data, and then takes the infinite volume limit when extracting the 
final results for the decay constants. We then take the larger of the following two values as 
our estimate of the finite volume error: 



VI. The difference between the central result and the result from a chiral fit in which the 
finite volume corrections are omitted. 



V2. The largest finite volume correction to the relevant data points, as determined by the 
central fit. For 0^)+ and 0b+; the "relevant data points" are the ones on each ensemble 
with the lightest valence mass, i.e., those closest to the chirally extrapolated point. 
For and 0^^) the relevant points are the ones on each ensemble with valence mass 
closest to m.. 



Method VI gives a larger difference for 0£)^ and 0b^; method V2 for 0^+ and 0^+ and the 
ratios. The resulting values are shown in Table XIII. Note that our choices are conservative 
because we correct for the (one-loop) finite volume errors, but nevertheless take the full size 
of these effects as our error. 
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VIII. RESULTS AND CONCLUSIONS 



After adding the error estimates described in the previous section in quadrature, we 
obtain: 



/^+ = 196.9(8.9) MeV, (8.1) 

fe, = 242.0(9.5) MeV, (8.2) 

IbJIb^ = 1.229(0.026), (8.3) 

fD+ = 218.9(11.3) MeV, (8.4) 

/d, = 260.1(10.8) MeV, (8.5) 

= 1.188(0.025). (8.6) 



Since our most rehable method of determining discretization errors combines them with 
statistical errors, we do not quote separate statistical and systematic errors. 

Figure 12 shows a comparison of our results for charmed decay constants with other lattice 
QCD calculations and with experiment. Our results agree with the only other three-flavor 
lattice QCD determination from the HPQCD collaboration [95], which is obtained with 
HISQ staggered valence quarks and asqtad staggered sea quarks. (The difference in f^^ is 
a bit greater than la.) They are also consistent with the two-flavor results of the ETM 
Collaboration using twisted-mass Wilson fermions [96], although the ETM error budget 
does not include an estimate of the uncertainty due to quenching the strange quark. One 
can also compare with "experimental" determinations of fo and if one assumes CKM 
unitarity to obtain the matrix elements \Vcd\ and |V^s|. For the D meson, Rosner and Stone 
combine CLEO's measurement of branching fraction B{D~^ — )■ /i"*"!^) [97] with the latest 
determination of \Vcd\ from the PDG [51] to obtain fo = 206.7(8.9) MeV [50]. For the 
meson, they average CLEO and Belle results for B{Df — )■ /U+i/) [98, 99] with CLEO and 
BABAR results for B{D^ — >■ t^i') [98, 100-102] to obtain a combined average for the two 
decay channels of fr,^ = 257.5(6.1) MeV [50]. The Heavy Flavor Averaging Group obtains 
a similar average, fc^ = 257.3(5.3) MeV [103]. Our results are consistent with these values, 
confirming Standard Model expectations at the ~ 5% level. 

Figure 13 shows a similar comparison of our results for bottom meson decay constants 
with other lattice QCD calculations. Our results agree with the published three-flavor 
determination using NRQCD 6-quarks and Asqtad staggered light quarks of the HPQCD 
collaboration [105], but are only marginally consistent with HPQCD's more recent calcu- 
lation of fs^ using HISQ light valence quarks [104]. Our results are also consistent with 
the two-flavor results of the ETM collaboration [96], who use Wilson heavy quarks and in- 
terpolate between the charm-mass region and the static limit to obtain results for bottom. 
Further, our result for the ratio fss/ Jb also agrees with the significantly less precise three- 
flavor determination using static 6-quarks and domain-wall light quarks by the RBC and 
UKQCD Collaborations [106]. 

For the D system the largest uncertainties in our current calculation stem from heavy- 
quark discretization, while the chiral extrapolation, the Zy factors, excited states, heavy- 
quark tuning, and the chiral-continuum extrapolation play important but subdominant roles. 
For the B system, heavy-quark tuning, statistics, and excited states are the sources of the 
largest errors, while the Zy factors and the chiral-continuum extrapolation (incorporating 
our estimate of heavy-quark discretization effects) are next in size. Recall that a novel 
feature of our work is the treatment of heavy-quark discretization effects, via the functions 
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fi in Eq. (2.4), and priors constraining the chiral-continuum fits to follow this form. At 
tree level, we have explicit calculations of the mismatch, some of which appeared already in 
Ref. [33] and all of which are compiled in Ref. [60]. Beyond the tree level, the continuum 
and static limits can be used to constrain the functional form. That said, the theoretical 
guidance of the priors cannot be highly effective in an analysis, such as this, with only two 
lattice spacings. Indeed, the quoted heavy-quark discretization errors are similar to less 
sophisticated power-counting estimates. 
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FIG. 12. Comparison of Jd and f^^ with other two- and three-flavor lattice QCD calculations and 
with experiment. Results shown come from Refs. [50, 95-102]. The HPQCD value is computed 
from their update to fo^ and their earlier result for the ratio fos/fo- 
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While completing this analysis, we have begun runs to generate data that will address 
the main sources of uncertainty reported here. The new data set will contain four times 
the configurations used here to reduce the statistical errors in the correlation functions and, 
thus, directly improve the decay amplitudes, the determinations of the hopping-parameters 
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FIG. 13. Comparison of Jb and Jb^ with other two- and three-flavor lattice QCD calculations. 
Results shown come from Refs. [96, 104-106]. In the case of Jb^ HPQCD has two separate calcu- 
lations using NRQCD b quarks and using HISQ b quarks; we show both the published NRQCD 
result (HPQCD '09) and the more recent HISQ result (HPQCD '11) in the plot above. 
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Kc and Kb, and the renormalization factors Zy^^ and ZyA^^ all of which feed into the decay 
constant. Our new data will also encompass two finer lattice spacings of a ~ 0.06 fm and 
a ~ 0.045 fm, in order to explicitly reduce light- and heavy-quark discretization errors and 
better control the continuum extrapolation. With four lattice spacings, our new method of 
heavy-quark discretization priors will be put to a more stringent test. The new runs will 
also include light valence- and sea-quark masses down to ~ ms/20 in order to better control 
the chiral extrapolation to the physical d and u quark masses. 

In order to reduce errors further, we will have to eliminate the errors from the matching 
factors and from quenching the charmed quark. The MILC Collaboration [107] is generating 
ensembles with 2+1+1 flavors of sea quarks with the HISQ action, with plans to provide 
a range of lattice spacings and sea quark masses equal to or more extensive than the 2+1 
asqtad ensembles. Use of the HISQ action for the charm valence quark will allow us to 
further reduce many of the uncertainties, and provides the particularly nice advantage that 
one can use the local pseudoscalar density without multiplicative renormalization to obtain 
the continuum matrix element [108]. In several years, once the full suite of HISQ ensembles 
with several sea-quark masses and lattice spacings has been analyzed, we expect to obtain 
percent-level errors for both B- and D-meson decay constants. This will enable precise tests 
of the Standard Model and may help to reveal the presence of new physics in the quark-flavor 
sector. 
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Appendix A: Heavy-quark Discretization Effects 

We are using the heavy-quark Lagrangian as given in [33] , with Kt = (or, equivalently, 
( = 1), = 1, and cb = ce = csw- This amounts to the Sheikholeslami-Wohlert La- 
grangian [32] for Wilson fermions [59]. The current has a heavy quark of this type, rotated 
as in Eq. (4.2) (c/. Eqs. (7.8)-(7.10) of Ref. [33]), and a staggered hght quark. At the tree 
level, the heavy-quark rotation is the same no matter what the other quark is. The dis- 
cretization effects are estimated from a (continuum) effective field theory [36-38], as shown 
explicitly for decay constants in Eqs. (8.7)-(8.12) of Ref. [36]. 

1. Theory 

Both QCD and lattice gauge theory can be described via 

i 

ClGT = >^^HQET(moa) = Cf'\mQ,moa)Oi, (A2) 

i 

where the Cj are short-distance coefficients and the Oi are operators describing the long- 
distance physics. The coefficients have dimension 4 — dimOj. For lattice gauge theory, 
they depend on niQa, which is a ratio of short distances a and l/mq. The effective-theory 
operators Oi in Eqs. (Al) and (A2) are the same. 
The error from each term is simply the difference 

error, = | [C^^H^Q, ^oa) - Cr*("^Q)] Oi\ . (A3) 

The relative error in our matrix elements can be estimated by setting {Oi) ~ Aq™j^'~^; 

choices for the QCD scale Aqcd are discussed below. The coefficient mismatch can be 
written 

C^\mQ, m,a) - Cr*("^Q) = a'^^''^-' h{moa). (A4) 

This recovers the usual counting of powers of a (familiar from Symanzik [34, 35]), but 
maintaining the full moa dependence. The final expression for the discretization errors is 
then 

error, oc U{m^a){aK^cBt^''^-\ (A5) 

For Wilson fermions, limmoa-s>o fi = constant (whereas in lattice NRQCD without fine tuning 
this is not the case). We have explicit calculations of the /, for the 0(a) and O(a^) errors at 
the tree level [33, 60]. The next subsection discusses how to use them to guide a continuum- 
limit extrapolation the 0{asa) and O(a^) errors. 

Equations (Al) and (A2) can be generalized to currents. For the axial- vector current, 

= C^^f{mQ)qzY^^,K - C^^^\mQ)v^q^,K -Y.Bm\^q)Q%. (A6) 

i 

Kt = C'2[imQ,moa)qtj^^-f5K - C'^^imQ,moa)v^'q-f5K - B'^liruQ , moa) Q''^,, (A7) 

i 

and = again means in the sense of matrix elements. Here selects the temporal component 
and ± the spatial, and the list of dimension-4 operators Q can be found in Refs. [37]. 
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The matrix element of the temporal component of the axial- vector current [cf. Eq. (4.7)] 
is normalized by multiplying with = C^°'^*/C]f*. The current mismatch then leads to 
errors 

^dim Q.-^f.^rnoa) = Z^.B'^l - B^, (A8) 
with the sum running over the two operators Q that point in the temporal direction [37]. 



2. Error Estimation 

The total error from heavy-quark discretization effects is then 

error = ^Zi (aAQCD)'7i("^oa) (A9) 

i 

where the sum runs over Lagrangian operators Oi of dimension 5 and 6 and current operators 
Qi of dimension 4 and 5, = dim (9j — 4 or dim Qi — 3, and the Zi are unknown coefficients. 
The functions fi (summarized below) have been computed for O(a^) and estimated for 
O(asa). We omit contributions of order a^a^, whether from extra operators or from iterating 
to second order operators with coefficients of order a^a. 

In the past, we have taken a very conservative Aqcd = 700 MeV and assumed a Gaussian 
distribution for the Zi centered on and of width 1. This amounts to treating the discretiza- 
tion errors as independent and adding them in quadrature. It also implicitly assumes that 
the data have nothing to say about the size or relative importance of the terms. 

Here, however, we incorporate these errors into the chiral-continuum extrapolation, dis- 
cussed in Sec. VI. This means that the Zi are now constrained fit parameters, with prior 
constraints discussed in Sec. VI. 

The fi are collected next. 



errors 



We start with these, because explicit expressions for the functions fi{mQa) are available. 
The Lagrangian leads to two bilinears, hD ■ Eh and /liS ■ [D x E]h, and many four-quark 
operators. At the tree level the coefficients of all four-quark operators vanish. At the tree 
level the coefficients of the two bilinears are the same, and the mismatch function is 



fE{moa) 



Using explicit expressions for l/m2 [33] and l/m\ [60], one finds 



We are using ce = 1, so 



1 



Ce{1 + niQa) - 1 
moa(2 + moa)(l + mod) 4(l + moa)" 



fE{moa) 



2 + SniQa 



8(2 + moa){l + moa)^ 



(AlO) 



(All) 



(A12) 



With no further assumptions, this term enters twice independently, so we take the width of 
this prior to be \/2 rather than 1. 
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The current leads to three more terms with non-zero coefficients, qTD'^h, gFiS • Bh, and 
qTa ■ Eh, which can be deduced from Eq. (A17) of Ref. [33]. Their coefficients can be read 
off from Eq. (A19). When cb = Ts the ffist two share the same coefficient 



/x(moa) 



(di{l + moo) 1 

moa(2 + moo) 2(2m2a)^' 

1 1 

+ 



(2 + moa)(l + moa) 
1 



2(1 + rriQa) 



4(l + moa)2 (2 + moa)2j' 



2(1 + moa) \2(2 + moa)(l + mofl) 



(A13) 



where the last term on the second line comes from using the tree-level di (as we do in the 
simulations). Because of the two-fold appearance, we again take the prior width to be a/2. 
For qTa ■ Eh 



/y(moa) 



di ({1 - ce){1 + moa) 



m2a moa{2 + moa) 
2 + 4moa + (mga)^ 



4(1 + moa)2(2 + moa)2' 
where the last line reflects the choices made for ce and di. 



(A14) 



b. 0(0^0) and O(a^) errors 

Here the mismatch functions fi{moa) start at order a^, and we do not have explicit ex- 
pressions for them. We take unimproved tree-level coefficients as a guide to the combinatoric 
factors and the asymptotic behavior as moa — )■ and moa — >■ 00. 

The Lagrangian leads to two bilinears, the kinetic energy O2 = hD'^h and the chromo- 
magnetic moment Ob = hiH ■ Bh. We match the former nonperturbatively, by identifying 
the meson's kinetic mass with the physical mass; the discretization error /2 stems, therefore, 
from discretization effects in M2. 

The computed kinetic meson mass is 



where [55] 



M2 = m2{K) + continuum binding energy + 5M2, (A15) 

A2 



Gmr 



m-i 



Tfl 

1 ) + Awi{m2ay 



(A16) 



and m2, ^4, and W4 are functions of moa and, hence, k. (See Refs. [33, 60] for explicit 
expressions.) Equating M2 to a physical meson mass means that we choose n such that 
m2{K) + SM2 = mq, thereby making in a relative error 

T M ^ \ J f 1 1 \ j6M2 , , 

error2 = A = A ^ ^ A — f . A17 

\2m2 2mQ) \2mQ-25M2 2mQ ) 2m| ^ ' 

The right-most expression is (aA)'^ /2(moa), /2 = [ ]/12(m2a)^, where [ ] is the bracket in 
Eq. (A16). It is formally smaller than the other errors considered here — /2 is of order 1 for 
all moa. Numerically, however, it is not much smaller. 
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At the tree level the chromomagnetic mismatch is 



fS\^oa) = -^ r. (A18) 

2(1 + rrioaj 

This has the right asymptotic behavior in both limits, so our Ansatz for the one-loop mis- 
match function is simply 

fBimoa) = ' r, (A19) 

2(1 -I- moo) 

and error^ is this function multiplied by aA. We take = 0.288 on the a ^ 0.12 fm 
ensembles, which is the value determined for ay from the plaquette [76] with one-loop 
running to scale q* = 2.5/a. On other ensembles, is found by assuming that the measured 
average taste splitting goes like a^a^ (with a determined from ri/a). This gives values 
that track av{q* = 2.5/a) quite well, which is why we make that q* choice. The results are 
rather insensitive to the details here. For example, using as = 0.325 on the a ^ 0.12 fm 
ensembles, which corresponds to q* = 2.0/a, increases the error estimate by less than 0.6 
MeV for /^+, and less than 0.25 MeV for 

The current leads to one more term, with tree-level mismatch function 



f? (moa) = ^ - rfi , ( A20) 

■'^ ^ ^ ^ 2(2 + moa)(l +moa) ^ ^ 

and the tree- level di is chosen so that = 0. As with the mismatch function Jb, we would 
like to anticipate fj^^ by setting S^'^ = and multiplying the rest with a^. But it is not 
generic that this vanishes as moct — ?• 0. Therefore, we take 

hM = -— ^ r, (A21) 

2(2 + mQa) 

which has the right asymptotic behavior. We take the prior width as \/2, because A'^ has 
two such corrections [37]. 



3. Dispersion relation, Eq. (3.1) 

We take a similar approach to the dispersion relation, Eq. (3.1), with the difference that 
we now know the sign of the leading effect. 
The tree-level functions are 

(A22) 
(A23) 

The binding energy enters and A^/ via the meson's kinetic energy. Hence, the binding 
contributions are 

A' = 2 1_ 

4 [0] , [0] / [0] NO / [0] no' K^^^J 

a (^4 aj'^ [1712 a) [m^a)-^ 
\ m i) a 



[0] 

ay 


1 


m[ 


''a 










(m2°'a)2 


I [0] 

[my 


'a) 


[0] 


[0] [0] 

= m\ awy . 
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and in Eq. (3.4) the binding energy floats within a Gaussian prior described by (A, 0"^) = 
(600,400) MeV. This choice conservatively brackets the binding energy of a heavy-strange 
meson. For the higher-order perturbative contribution to the coeflicients, we take the 
Ansatze based on the asymptotic behavior: 



2/4 + -24 ln(l -|- mga) 



1 + moaY 
[1] Vi'i^Q^L + z^i ln(l + rriQa) 
1 + rriQa 



(A26) 



0^4' — ^ , _ ; l^^'J 



where the ys and 2;s float within Gaussian priors described by {y4,ay^) = (3,5), {z4^,az4^) = 
(1,2), {y4',<Jy^,) = (0,0), and {z4',az^,) = (0,2). The terms proportional to yi stem from 
the rriQa — )■ limit, in which the renormalization of must coincide with that of mi, and 
^4 = iTiiaWi must vanish like moa. The terms proportional to Zi stem from the niQa — )■ oo 
limit, where the static limit is obtained. Except for 1/4/, the numerical values have been 
chosen consistent with one-loop experience for nii and m2 [109]. We have set 1/4/ = 0, 
because at small m^a it is indistinguishable from the other term in a^,', and our range of 
rriQa does not reach far into the region moct ^ 1. 
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Appendix B: Two point fit results from Analysis I 



TABLE XIV. Heavy-light pseudoscalar meson masses and renormalized decay amplitudes obtained 
from Analysis I fits of the charm correlators at lattice spacing a w 0.09 fm. 



ami jams CLf^q aMn cC'I'^^h x^/dof p 






0031/0 


031 





0031 





7523 


;o 


0016) 





0857 


;o 


0015) 


58/48 





23 





0031/0 


031 





0044 





7553 


;o 


0014) 





0873 


;o 


0013) 


56/48 





28 





0031/0 


031 





0062 





7589 


;o 


0011) 





0890 


;o 


0011) 


55/48 





33 





0031/0 


031 





0087 





7634 


;o 


0009) 





0910 


;o 


0009) 


53/48 





38 





0031/0 


031 





0124 





7699 


;o 


0007) 





0936 


;o 


0007) 


53/48 





41 





0031/0 


031 





0186 





7807 


;o 


0005) 





0978 


;o 


0006) 


52/48 





44 





0031/0 


031 





0272 





7954 


;o 


0004) 





1030 


;o 


0005) 


50/48 





5 





0031/0 


031 





031 





8018 


;o 


0004) 





1052 


;o 


0004) 


50/48 





5 





0062/0 


031 





0031 





7541 


;o 


0030) 





0875 


;o 


0027) 


56/48 





37 





0062/0 


031 





0044 





7577 


;o 


0023) 





0899 


;o 


0021) 


52/48 





49 





0062/0 


031 





0062 





7613 




0019) 





0917 




0018) 


50/48 





58 





0062/0 


031 





0087 





7654 




0015) 





0933 




0015) 


58/51 





43 





0062/0 


031 





0124 





7712 


;o 


0012) 





0952 


;o 


0012) 


52/48 





48 





0062/0 


031 





0186 





7810 


;o 


0009) 





0985 


^0 


0010) 


56/48 





37 





0062/0 


031 





0272 





7952 


;o 


0006) 





1032 


;o 


0008) 


59/48 





28 





0062/0 


031 





031 





8015 


;o 


0005) 





1052 


;o 


0007) 


60/48 





25 





0124/0 


031 





0031 





7551 


;o 


0038) 





0930 


;o 


0036) 


60/48 





27 





0124/0 


031 





0042 





7554 


;o 


0031) 





0926 


;o 


0028) 


65/48 





15 





0124/0 


031 





0062 





7574 




0023) 





0929 


;o 


0021) 


65/48 





16 





0124/0 


031 





0087 





7608 


;o 


0017) 





0938 


;o 


0015) 


59/48 





28 





0124/0 


031 





0124 





7666 


;o 


0013) 





0957 


;o 


0012) 


49/48 





63 





0124/0 


031 





0186 





7766 


;o 


0008) 





0991 


;o 


0009) 


42/48 





85 





0124/0 


031 





0272 





7907 


;o 


0006) 





1038 


;o 


0007) 


48/48 





64 





012-1/0 


031 





031 





79G9 


;o 


0005) 





1058 


;o 


OOOG) 


53/18 





47 
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TABLE XV. Heavy-light pseudoscalar meson masses and renormalized decay amplitudes obtained 
from Analysis I fits of the charm correlators at lattice spacing a 0.12 fm. 



ami/ am s o.'f^q ciMh a^/^(j)H x^/dof p 






005/0 


050 





005 





9943 


;o 


0032) 





1436 


;o 


0030) 


30/30 





52 





005/0 


050 





007 





9977 


;o 


0024) 





1453 


;o 


0024) 


29/30 





6 





005/0 


050 





01 


1 


0026 




0018) 





1477 


'0 


0019) 


28/30 





64 





005/0 


050 





014 


1 


0090 


;o 


0016) 





1508 


^0 


0017) 


28/30 





64 





005/0 


050 





02 


1 


0186 


;o 


0013) 





1551 


;o 


0015) 


29/30 





58 





005/0 


050 





03 


1 


0345 


;o 


0010) 





1620 


;o 


0012) 


33/30 





42 





005/0 


050 





0415 


1 


0526 


;o 


0008) 





1694 


;o 


0010) 


36/30 





27 





007/0 


050 





005 





9948 


;o 


0035) 





1442 


;o 


0035) 


17/30 





98 





007/0 


050 





007 





9975 




0027) 





1455 


'0 


0028) 


19/30 





95 





007/0 


050 





01 


1 


0019 


;o 


0021) 





1476 


^0 


0021) 


22/30 





89 





007/0 


050 





014 


1 


0081 


;o 


0016) 





1504 


^0 


0017) 


24/30 





83 





007/0 


050 





02 


1 


0178 


;o 


0012) 





1547 


;o 


0014) 


23/30 





85 





007/0 


050 





03 


1 


0338 


;o 


0009) 





1615 


;o 


0010) 


20/30 





94 





007/0 


050 





0415 


1 


0520 


;o 


0007) 





1687 


;o 


0008) 


19/30 





95 





010/0 


050 





005 





9958 


;o 


0039) 





1461 


;o 


0041) 


15/30 





99 





010/0 


050 





007 


1 


0000 




0031) 





1486 


^0 


0032) 


20/30 





94 





010/0 


050 





01 


1 


0057 


;o 


0024) 





1516 


^0 


0026) 


26/30 





75 





010/0 


050 





014 


1 


0126 


;o 


0019) 





1549 


^0 


0021) 


29/27 





41 





010/0 


050 





02 


1 


0226 


;o 


0015) 





1594 


;o 


0017) 


33/30 





39 





010/0 


050 





03 


1 


0387 


;o 


0011) 





1662 


;o 


0014) 


31/30 





5 





010/0 


050 





0415 


1 


0567 


;o 


0008) 





1733 


;o 


0011) 


27/30 





68 





020/0 


050 





005 





9942 


;o 


0046) 





1537 


;o 


0050) 


49/30 





036 





020/0 


050 





007 





9959 




0036) 





1533 


^0 


0039) 


49/30 





036 





020/0 


050 





01 





9987 


;o 


0027) 





1532 


^0 


0031) 


48/30 





051 





020/0 


050 





014 


1 


0037 


;o 


0021) 





1543 


^0 


0024) 


45/30 





075 





020/0 


050 





02 


1 


0124 


;o 


0016) 





1575 


;o 


0019) 


43/30 





11 





020/0 


050 





03 


1 


0274 


;o 


0011) 





1632 


;o 


0014) 


37/30 





27 





020/0 


050 





0415 


1 


0447 


;o 


0009) 





1695 


;o 


0012) 


32/30 





48 





030/0 


050 





005 





9830 


;o 


0042) 





1475 


;o 


0042) 


33/30 





39 





030/0 


050 





007 





9853 




0033) 





1485 


'0 


0033) 


33/30 





4 





030/0 


050 





01 





9897 


;o 


0025) 





1505 


^0 


0025) 


32/30 





47 





030/0 


050 





014 





9960 


;o 


0020) 





1534 


^0 


0020) 


31/30 





53 





030/0 


050 





02 


1 


0054 


;o 


0015) 





1574 


;o 


0016) 


32/30 





46 





030/0 


050 





03 


1 


0205 


;o 


0011) 





1633 


;o 


0012) 


37/30 





27 





030/0 


050 





0415 


1 


0376 


;o 


0009) 





1695 


;o 


0010) 


40/30 





15 
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TABLE XVI. Heavy-light pseudoscalar meson masses and renormalized decay amplitudes obtained 
from Analysis I fits of the charm correlators at lattice spacing a 0.15 fm. 



ami 1 o-iT^s 




amq 


aMu 








xVdof p 




0.0097/0 


0484 


0.0048 


1.1659(0 


0044) 


0.1979(0 


0052) 


20/20 





5 


0.0097/0 


0484 


0.007 


1.1710(0 


0034) 


0.2017(0 


0040) 


22/20 





37 


0.0097/0 


0484 


0.0097 


1.1768(0 


0027) 


0.2054(0 


0032) 


25/20 





26 


0.0097/0 


0484 


0.0194 


1.1951(0 


0016) 


0.2159(0 


0020) 


25/20 





26 


0.0097/0 


0484 


0.029 


1.2117(0 


0012) 


0.2242(0 


0015) 


20/20 





51 


0.0097/0 


0484 


0.0484 


1.2432(0 


0009) 


0.2385(0 


0012) 


15/20 





79 


0.0194/0 


0484 


0.0048 


1.1726(0 


0046) 


0.2106(0 


0052) 


23/20 





35 


0.0194/0 


0484 


0.007 


1.1749(0 


0036) 


0.2105(0 


0041) 


23/20 





35 


0.0194/0 


0484 


0.0097 


1.1785(0 


0028) 


0.2113(0 


0031) 


23/20 





32 


0.0194/0 


0484 


0.0194 


1.1935(0 


0016) 


0.2174(0 


0020) 


30/20 





092 


0.0194/0 


0484 


0.029 


1.2091(0 


0013) 


0.2244(0 


0016) 


32/20 





055 


0.0194/0 


0484 


0.0484 


1.2400(0 


0010) 


0.2381(0 


0013) 


27/20 





17 


0.0290/0 


0484 


0.0048 


1.1613(0 


0044) 


0.1975(0 


0049) 


17/20 





72 


0.0290/0 


0484 


0.007 


1.1660(0 


0034) 


0.2010(0 


0040) 


18/20 





64 


0.0290/0 


0484 


0.0097 


1.1717(0 


0026) 


0.2049(0 


0031) 


21/20 





47 


0.0290/0 


0484 


0.0194 


1.1896(0 


0015) 


0.2151(0 


0019) 


24/20 





3 


0.0290/0 


0484 


0.029 


1.2058(0 


0011) 


0.2229(0 


0015) 


23/20 





32 


0.0290/0 


0484 


0.0484 


1.2368(0 


0008) 


0.2364(0 


0011) 


20/20 





49 
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TABLE XVII. Heavy- light pseudoscalar meson masses and renormalized decay amplitudes obtained 
from Analysis I fits of the bottom correlators at lattice spacing a ^ 0.09 fm. 



ami/ am s aMjj a^/^^ij x^/dof p 






0031/0 


031 





0031 


1 


6509 


;o 


0018) 





1359 


;o 


0016) 


41/39 





48 





0031/0 


031 





0044 


1 


6532 


;o 


0016) 





1378 


;o 


0015) 


40/39 





51 





0031/0 


031 





0062 


1 


6562 




0015) 





1402 




0014) 


40/39 





49 





0031/0 


031 





0087 


1 


6601 


;o 


0013) 





1433 


;o 


0014) 


42/39 





42 





0031/0 


031 





0124 


1 


6659 


;o 


0012) 





1475 


;o 


0013) 


45/39 





31 





0031/0 


031 





0186 


1 


6752 


;o 


0011) 





1542 


;o 


0012) 


47/39 





23 





0031/0 


031 





0272 


1 


6879 


;o 


0009) 





1628 


;o 


0011) 


49/39 





19 





0031/0 


031 





031 


1 


6934 




0009) 





1664 




0011) 


49/39 





18 





0062/0 


031 





0031 


1 


6539 


;o 


0046) 





1358 


;o 


0051) 


40/39 





56 





0062/0 


031 





0044 


1 


6557 


;o 


0039) 





1377 


;o 


0044) 


37/39 





68 





0062/0 


031 





0062 


1 


6584 


;o 


0032) 





1402 


;o 


0037) 


34/39 





77 





0062/0 


031 





0087 


1 


6620 


;o 


0027) 





1434 


;o 


0031) 


34/39 





8 





0062/0 


031 





0124 


1 


6675 


;o 


0022) 





1480 


;o 


0026) 


36/39 





72 





0062/0 


031 





0186 


1 


6767 


;o 


0018) 





1550 


;o 


0022) 


41/39 





53 





0062/0 


031 





0272 


1 


6892 


;o 


0014) 





1637 


;o 


0019) 


45/39 





37 





0062/0 


031 





031 


1 


6946 


;o 


0014) 





1672 


;o 


0018) 


45/39 





35 





0124/0 


031 





0031 


1 


6532 




0036) 





1387 


;o 


0038) 


52/39 





16 





0124/0 


031 





0042 


1 


6550 


;o 


0033) 





1407 


;o 


0034) 


48/39 





27 





0124/0 


031 





0062 


1 


6576 


;o 


0030) 





1432 


;o 


0031) 


40/39 





55 





0124/0 


031 





0087 


1 


6606 


;o 


0027) 





1456 


;o 


0029) 


35/39 





77 





0124/0 


031 





0124 


1 


6650 


;o 


0024) 





1488 


;o 


0027) 


33/39 





84 





0124/0 


031 





0186 


1 


6730 


;o 


0019) 





1544 




0023) 


36/39 





73 





0124/0 


031 





0272 


1 


6847 


;o 


0016) 





1623 




0021) 


42/39 





48 





0124/0 


031 





031 


1 


6900 


;o 


0015) 





1657 


;o 


0020) 


45/39 





38 



57 



TABLE XVIII. Heavy-light pseudoscalar meson masses and renormalized decay amplitudes ob- 
tained from Analysis I fits of the bottom correlators at lattice spacing a 0.12 fm. 



ami/ am s o.i^q aMn a^/^(f>H x^/dof p 






005/0 


050 





005 


1 


9170 


;o 


0044) 





2236 


'0 


0050) 


45/27 





03 





005/0 


050 





007 


1 


9197 


;o 


0039) 





2263 


;o 


0046) 


46/27 





022 





005/0 


050 





01 


1 


9235 


'0 


0033) 





2300 




0040) 


46/27 





021 





005/0 


050 





014 


1 


9287 


^0 


0029) 





2347 




0036) 


45/27 





027 





005/0 


050 





02 


1 


9367 


;o 


0024) 





2418 


;o 


0031) 


43/27 





046 





005/0 


050 





03 


1 


9503 


;o 


0020) 





2532 




0026) 


39/27 





096 





005/0 


050 





0415 


1 


9657 


;o 


0017) 





2654 


;o 


0023) 


36/27 





17 





007/0 


050 





005 


1 


9147 


;o 


0036) 





2224 




0039) 


37/27 





12 





007/0 


050 





007 


1 


9177 


^0 


0033) 





2254 




0037) 


35/27 





17 





007/0 


050 





01 


1 


9219 


^0 


0030) 





2292 


;o 


0036) 


34/27 





2 





007/0 


050 





014 


1 


9272 


^0 


0028) 





2337 


;o 


0037) 


35/27 





19 





007/0 


050 





02 


1 


9351 


;o 


0026) 





2401 


;o 


0037) 


36/27 





15 





007/0 


050 





03 


1 


9485 


;o 


0022) 





2508 


;o 


0035) 


38/27 





096 





007/0 


050 





0415 


1 


9638 


;o 


0019) 





2628 




0031) 


40/27 





07 





010/0 


050 





005 


1 


9182 


;o 


0047) 





2254 


;o 


0047) 


30/27 





4 





010/0 


050 





007 


1 


9207 


^0 


0041) 





2284 




0042) 


32/27 





29 





010/0 


050 





01 


1 


9250 


'0 


0035) 





2328 


;o 


0037) 


36/27 





18 





010/0 


050 





014 


1 


9307 


^0 


0030) 





2383 


;o 


0033) 


39/27 





097 





010/0 


050 





02 


1 


9391 


;o 


0025) 





2457 


;o 


0028) 


43/27 





048 





010/0 


050 





03 


1 


9527 


;o 


0020) 





2569 


;o 


0024) 


47/27 





02 





010/0 


050 





0415 


1 


9682 


;o 


0017) 





2689 


;o 


0021) 


51/27 





0092 





020/0 


050 





005 


1 


9136 


;o 


0060) 





2278 


;o 


0069) 


33/27 





27 





020/0 


050 





007 


1 


9163 


^0 


0050) 





2305 




0059) 


33/27 





28 





020/0 


050 





01 


1 


9200 


^0 


0042) 





2340 


;o 


0050) 


31/27 





36 





020/0 


050 





014 


1 


9249 


^0 


0036) 





2381 


;o 


0043) 


29/27 





47 





020/0 


050 
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TABLE XIX. Heavy-light pseudoscalar meson masses and renormalized decay amplitudes obtained 
from Analysis I fits of the bottom correlators at lattice spacing a 0.15 fm. 
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